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Abstract 

The  distribution  iteration  (DI)  algorithm,  developed  by  Wager  [32]  and  Prins 
[28],  for  solving  the  Boltzmann  Transport  Equation  (BTE)  has  proven,  with  further 
development,  to  be  a  robust  alternative  to  von  Neumann  iteration  on  the  scattering 
source,  aka  source  iteration  (SI).  Previous  work  with  DI  was  based  on  the  time- 
independent  form  of  the  transport  equation.  In  this  research,  the  DI  algorithm  was 

1.  Improved  to  provide  faster,  more  efficient,  robust  convergence 

2.  Extended  to  XYZ  geometry 

3.  Extended  to  Multigroup  Energy  treatment 

4.  Extended  to  solve  the  time-dependent  form  of  the  Boltzmann  Transport  Equa¬ 
tion. 

The  discrete  ordinates  equations  for  approximating  the  BTE  have  been  solved 
using  SI  since  the  discrete  ordinates  method  was  developed  at  Los  Alamos  Scientific 
Laboratory  by  1953.  However,  SI  is  often  inefficient  by  itself  and  requires  an  accelera¬ 
tor  in  order  to  produce  results  efficiently  and  reliably.  The  acceleration  schemes  that 
are  in  use  in  production  codes  are  Diffusion  Synthetic  Acceleration  (DSA)  and  Trans¬ 
port  Synthetic  Acceleration  (TSA).  DSA  is  ineffective  for  some  problems,  and  cannot 
be  extended  to  high-performance  spatial  quadratures.  TSA  is  less  effective  than  DSA 
and  fails  for  some  problems.  Krylov  acceleration  has  been  explored  in  recent  years, 
but  has  many  parameters  that  require  problem-dependent  tuning  for  efficiency  and 
effectiveness. 

The  DI  algorithm  is  an  alternative  to  source  iteration  that,  in  our  testing,  does 
not  require  an  accelerator.  I  developed  a  formal  verification  plan  and  executed  it  to 
verify  the  results  produced  by  my  code  that  implemented  DI  with  the  above  features. 
A  new,  matrix  albedo,  boundary  condition  treatment  was  developed  and  implemented 


IV 


so  that  infinite-medium  benchmarks  could  be  included  in  the  verification  test  suite. 
The  DI  algorithm  was  modified  for  parallel  efficiency  and  the  prior  instability  of  the 
refinement  sweep  was  corrected.  The  testing  revealed  that  DI  performed  as  well  or 
faster  than  source  iteration  with  DSA  and  that  DI  continued  to  work  where  DSA 
failed.  Performance  did  degrade  when  the  diamond-difference  (without  fixup)  spatial 
quadrature  was  used.  Because  diamond-difference  is  a  non-positive  spatial  quadra¬ 
ture,  it  can  produce  nonphysical  negative  fluxes,  particularly  in  higher  dimensions. 
I  developed  a  new  fixup  scheme  to  accommodate  the  negative  fluxes,  but  it  did  not 
improve  performance  in  XYZ  geometry  when  the  scattering  ratio  was  near  unity. 

My  DI  algorithm  successfully  solves  the  time-dependent  form  of  the  BTE  using 
the  semi-implicit  method  implemented  by  PARTISN.  The  agreement  between  DI  and 
PARTISN  was  excellent. 

With  these  improvements  and  tests,  DI  is  ready  for  use  as  a  general  replacement 
for  von  Neumann  iteration  on  the  scattering  source. 
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Time  Dependent  Discrete  Ordinates  Neutron  Transport 
Using  Distribution  Iteration 
in  XYZ  Geometry 


I.  Introduction 

The  linearized  form  of  the  Boltzmann  Transport  Equation  (BTE),  aka  the  neu¬ 
tron  transport  equation,  can  be  expressed  as 


~ld_ 
v  dt 


+  O  •  V  +  cr(r,  E ,  t) 


ili(r,Ct,E,t)  =  q(f,(t,E,t), 


(1.1) 


where  r  is  the  position;  O  is  the  direction  of  travel;  E  is  the  energy  of  the  neutron; 
t  is  time;  v  is  the  neutron  velocity;  a  is  the  total  cross  section  and  ij}  is  the  angular 
flux.  The  rate  density,  q,  can  be  represented  as 


q(f,Cl,E,t)=  /  dE'  /  dU  crs(f,  E'  — >  E,  Cl'  ■  O,  t)i()(f,  Cl1,  E' ,  t) 


+  X(E)  /  d E' uaf{r,E',t)  /  dU  ^(r,  Cl1,  E\  t)  +  ge(r,  Cl,  E,  t),  (1.2) 


where  as  is  the  differential  scattering  cross  section;  y  is  the  fission  spectrum;  v  is  the 
mean  number  of  neutrons  produced  per  fission;  <jf  is  the  fission  cross  section  (Lewis 
&  Miller  [22])  and  qe  is  the  fixed  (embedded)  emission  rate  density.  The  rate  density 
can  be  further  categorized  as 

/»oo 

,E,t)=  d  E' 

Jo 

qf(r,£l,E,t)  =  x(E)  J  d E'  v<jf(r,  E' ,t)  j  dD'  ip(f,  Clf,  E',  t)  (1.4) 

where  qs  is  the  scattering  source  and  qJ  is  the  fission  source. 

The  vast  majority  of  all  research  in  neutron  transport  is  focused  either  on  the 
time-independent  BTE  or  on  Monte  Carlo  methods.  The  time-independent  BTE 


qs(r,  D 


J  dU  crs(f,  E'  — *  E,  Cl'  ■  Cl,  t)ip(r,  Cl',  E' ,  t)  (1.3) 


1 


simplifies  the  transport  equation  by  solving  for  the  steady-state  solution,  which  forces 
the  time  derivative  to  zero.  A  common  approach  for  solving  the  time-dependent  BTE 
utilizes  a  method  based  on  solving  the  time- independent  BTE  (Dupree  [11],  Reed  [30] 
and  Hill  [16]).  Monte  Carlo  methods  are  also  used  to  solve  time-dependent  problems; 
however,  they  can  require  a  large  number  of  particle  histories  in  order  to  reduce 
uncertainty. 

The  focus  of  this  research  was  to  incorporate  the  distribution  iteration  method 
into  a  multigroup,  time-dependent  BTE  solver.  The  distribution  iteration  algorithm 
has  substantial  advantages  over  the  more  common  source  iteration  algorithm,  the 
most  notable  ones  being  its  rapid  convergence  in  highly-scattering  regions  and  the 
ability  to  easily  implement  distribution  iteration  in  parallel. 

1.1  Background 

The  integro-differential  form  of  the  BTE  makes  it  particularly  difficult  to  solve 
and  analytic  solutions  are  only  possible  for  the  simplest  of  problems.  Numerical 
solutions  to  the  BTE  are  obtained  by  applying  a  variety  of  approximations  in  order 
to  simplify  the  BTE  (Carlson  [7]): 

•  Assume  a  stationary  (time-independent)  solution 

•  Discretize  in  energy  (multigroup  approximation) 

•  Discretize  in  angle  (discrete  ordinates) 

•  Spherical  harmonics  expansion  of  the  scattering  source 

The  direct  application  of  difference  methods  directly  to  the  BTE  is  usually  avoided 
because  neutrons  are  not  strictly  conserved  (Carlson  [9]).  In  order  to  ensure  neu¬ 
tron  conservation,  cell  averages  of  the  distribution  function  are  used  and  cell  balance 
equations  are  coupled  to  the  system  of  equations. 

The  spherical  harmonics  expansion  of  the  scattering  source  is  the  result  of  using 
a  Legendre  expansion  of  the  differential  cross  section  as  in  the  scattering  source  qs 


2 


making  the  expansion 


a9(r,  E'  ^E,Cl'  ■  Cl,  t )  =  ^(21  +  1  )asl(r,  E'  ->  E,  t)Pt(Cl'  ■  Cl),  (1.5) 

1=0 

where  P;  is  the  Legendre  polynomial  of  order  l  and  asi  are  the  Legendre  scattering 
moments.  Substituting  this  approximation  into  the  scattering  source  yields, 


q*(r,  Cl,  E,  t )  =  '$2(21  +  1)  /  d E'  asl(r ,  E'  E,  t ) 


1=0 


d Q!  Pl(Cl-Cl,)il)(r,Q!,E,,t).  (1.6) 


where  the  Legendre  expansion  has  been  truncated  to  L  terms.  Using  the  Legendre 
addition  theorem, 

1  1 

m  ■  ft)  =  E  itafnirtatsV),  (i.r) 

m=—l 

where  Y)m  is  the  spherical  harmonics  function  and  Yt*n  is  its  complex  conjugate,  per¬ 
mits  the  removal  of  the  inconvenient  scattering  angle  expression,  Cl'  ■  Cl,  and  replaces 
it  with  the  more  convenient  directions  of  travel  Cl'  and  Cl.  With  the  above  derivations, 
the  scattering  term  can  be  expressed  as 

l  i  r 

qs(r,Cl,E,t)  =  /  dE'  Esi(r,E'  -  E,t) 

1=0  m=—l  ' 

J  dn'Y^Ci'^CiCECt).  (i.8) 

The  angular  flux  coefficients  are  dehned  as 

4im(r,E',t)  =  J  dfi'  Ylm(Q')^(r,Cl',E',t),  (1.9) 
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which  yields 


l  i  r 

qs(  f,Cl,E,t)  =  J2Y1  /  d E'  cr si  if  \  E'  -  E,  t)(f>lm(r,  E',  t).  (1.10) 

1=0  m=—l  ^ 

This  form  of  the  scattering  term  when  combined  with  an  angular  quadrature,  the 
discrete  ordinates  approximation,  results  in  a  finite  set  of  equations  where  the  an¬ 
gular  dependence  has  been  explicitly  eliminated.  Given  an  angular  quadrature  of  N 
ordinates,  the  angular  flux  coefficients  can  be  determined  from 


Vh  (r,E',t) 
il> N(r,E',t ) 


(1.11) 


where  the  matrix  on  the  left  has  rows  of  l  and  m  and  Wn  are  the  angular  quadrature 
weights.  In  slab  geometry,  the  spherical  harmonics  functions  reduce  to  Legendre 
polynomials  clue  to  symmetry.  In  XY  geometry,  the  spherical  harmonics  functions 
reduce  to  just  the  even  terms.  However,  in  XYZ  geometry  the  spherical  harmonic 
functions  do  not  simplify. 

Having  discretized  the  angular  dependence,  we  now  seek  to  discretize  the  en¬ 
ergy  dependence.  The  energy  discretization  approach  used  in  almost  all  deterministic 
methods  is  the  multigroup  approximation,  which  divides  the  continuous  energy  range 
into  a  finite  set  of  groups.  This  transforms  the  integration  over  energy  into  a  summa¬ 
tion  over  the  energy  groups.  The  multigroup  approximation  defines  the  group  angular 
flux  ('ijjg)  as 

ipg(r,£l,t)  —  f  dE  ip(r,  Cl,  E,  t)  —  f  cl E  E,t),  (1.12) 

J  Eg  J  g 

where  g  is  the  energy  group  index  such  that  g  E  [1,  G],  where  G  is  the  total  number 
of  energy  groups.  The  above  definition  transforms  the  BTE  into  a  coupled  set  of 
monoenergetic  transport  equations,  which  are  then  solved  by  iteration. 
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The  first  step  in  deriving  the  within-group  transport  equation  is  the  integration 
of  equation  (1.1)  over  an  energy  group,  g,  yielding 


-4  +  6-v 

vn  at 


L  u9 


d E  ip(f,  11,  E,  t)  +  /  d E  a(r,  E,  t)^(r,  12,  E,  t) 

i  Jg 

[  d E  f  d E'  f  dll'  <7s(f,  £'  ->  E,  Cl'  ■  Cl,  t)if>(r,  Cl',  E’,  t ) 

„/=1  Jg  J  g1  J 


G 


+  /  dT?  x(E)  /  d E'  vcrf(r,  E' ,t)  /  dll7  ip(r,  Cl',  E' ,  t) 

J g  g'=l'9  ^ 


d E  qe(r,Cl,E,t),  (1.13) 


where  vg  is  the  group  velocity.  Because  the  neutron  velocity,  v,  is  a  continuous 
function  of  energy  and  the  angular  flux  if)  is  integrable  in  the  same  domain,  then  the 
first  mean  value  theorem  for  integration  allows 


d  E 


1  <9-0 
v{E)  dt 


1 


(1.14) 


where  £  is  some  energy  within  the  group.  Group  cross-sections,  fission  spectrum,  and 
fixed  source  can  be  defined  as  follows  (Lewis  &  Miller  [22])  to  simplify  the  represen¬ 
tation 


ag(f,Cl,t) 


Xg 


uafg  ('G  t) 


Jg  d E  cr(r,  E,  Cl,  E,  t ) 

^ g(r,Cl,t ) 

[  dEX(E), 

Jg 

Ig  dE  vof{r,  E,  t )  /  d!2  #r,  Cl,  E,  t ) 
J  dll  'ipg(r,  Cl',  t ) 


(1.15) 

(1.16) 

(1.17) 


and 

Qg(r,Cl,t)=  f  dE  qe(r,Cl,E,t).  (1.18) 

Jg 

The  multigroup  group-to-group  scattering  cross  section,  as,  is  considerably  more  com¬ 
plicated  to  define.  The  typical  approach  is  to  apply  the  discrete  ordinates  approxi- 
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mation,  which  yields 


vigg'(r,t) 


fg  d E  fg,  d E'  asi(f,  E’  ->  E,  t)(j>im(r,  E' ,  t ) 

(? 


(1.19) 


where 

4>im(r ,  E’,  t)  =  J  <m  Ylm(Cl)^(f,  Cl,  Ef,  t )  (1.20) 

and 

4>img'(r,t )  =  J  dfi  Yim(Ci)ipgi(r,Cl,t).  (1.21) 

While  an  angular  dependence  to  the  multigroup  total  cross  section,  ag,  has  been 
introduced,  it  can  be  eliminated  by  using  a  spherical  harmonics  expansion.  Applying 
the  discrete  ordinates  approximation,  the  multigroup  form  of  the  BTE  is 


~  +  Cln-V  +  ag(f,t) 

L  V9  °l 


i/>g(r,Cln,t ) 


oo  l  G 

E  E  VJn»)E  G Igg'iT 

1=0  m=—l  g'= 1 

G 

+  Xg^2  t )  +  qeg(f,  Cln ,  t),  (1.22) 

g’= i 


where  <pg  is  the  group  scalar  flux.  Segregating  the  within-group  g  from  the  scattering 
source  yields  the  within-group  BTE1 


'^d_ 

Vg  Dt 


+  Cln  •  V  +  ag(f,  t ) 


^Pg(j" i  Clni  t'} 

L 

Ypl  +  1  )Pi(Cln)aigg(f,  t)(j)ig{r ,  f)  +  g9(r,  f), 

1=0 


(1.23) 


where  represents  all  the  contributions  to  group  g  clue  to  scatter  from  other  energy 
groups  (downscatter  and  upscatter),  fission,  and  the  fixed  emitter. 

1  The  within-group  fission  source  is  typically  not  segregated  and  is  treated  as  a  known  value 
which  is  re-evaluated  after  the  within-group  solution  is  determined  (Carlson  &  Lathrop  [9]). 
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1.1.1  Source  Iteration.  The  standard  method  for  solving  the  within-group 
BTE  is  the  scattering  source  iteration  (SI),  which  is  based  on  the  von  Neumann  series 
solution.  The  angular  flux,  ip,  can  be  decomposed  into  the  sum  of  collided  angular  flux 
contributions,  ip171.  The  'ipm  take  on  the  physical  significance  that  ip°  is  the  uncollided 
flux,  ip1  is  the  one-collision  flux,  etc.  We  can  define  a  series  of  iterates  (Lewis  & 
Miller  [22]),  such  that 

ip°  =  0 
ip1  =  4>° 

i 

=  (1.24) 

m= 0 

Thus,  ipl+1  is  the  flux  of  neutrons  that  have  had  l  or  fewer  collisions.  A  consequence 
of  this  method  is  that  a  problem  that  has  a  scattering-to-total  cross  section  ratio,  c, 
approaching  unity  requires  an  increasing  number  of  iterations  in  order  to  converge  to 
a  solution.  Furthermore,  while  the  individual  contribution  to  the  angular  flux  from  a 
specific  iterate  ipm)  where  m  is  large,  might  be  small,  the  sum  of  the  contributions  from 
the  high  scatter  iterates  is  not  necessarily  small.  This  can  cause  the  source  iteration 
method  to  falsely  converge  to  a  solution  by  terminating  the  iteration  prematurely. 
In  order  to  accelerate  convergence  to  a  solution,  a  variety  of  techniques  are  used, 
the  most  common  one  being  diffusion  synthetic  acceleration  (Alcouffe  [3]).  Using 
diffusion  synthetic  acceleration  with  source  iteration  is  a  well  established  method 
used  extensively  in  production  codes. 

Diffusion  synthetic  acceleration  generally  performs  well  for  a  broad  range  of 
problems;  however,  it  does  lose  effectiveness  if  the  problem  does  not  exhibit  sufficient 
diffusive  behavior.  Furthermore,  it  is  particularly  difficult  to  implement  the  newer, 
better  performing  spatial  quadratures  that  have  been  recently  developed  into  diffusion 
synthetic  acceleration.  While  other  accelerators  exist,  none  of  them  are  a  panacea. 
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1.1.2  Distribution  Iteration.  Distribution  iteration  (DI)  offers  several  ad¬ 
vantages  over  the  more  traditional  source  iteration  method,  most  notably  its  rapid 
convergence  when  the  scattering-to-total  ratio  approaches  unity  and  its  inherently 
parallel  nature  (Wager  [32]  and  Prins  [28]).  The  central  idea  behind  DI  is  to  reduce 
the  global  transport  problem  into  a  coupled-cell  partial  current  problem  that  can  be 
solved  directly.  The  basic  DI  algorithm  is  shown  in  Algorithm  1.1.  DI  differs  from 
source  iteration  in  that  it  solves  the  contribution  due  to  scattering  directly  and  does 
not  iterate  on  the  scattering  source.  This  key  difference  eliminates  the  need  to  use  an 
accelerator  to  reduce  the  number  of  iterations  required  for  convergence. 

Guess  the  angular  flux  distribution  on  each  face,  e.g.  isotropic 
Compute  the  within  group  source  due  to  fixed  emission  sources 

repeat 

Use  the  angular  flux  distributions  to  determine  the  face-to-face  transport  prob¬ 
abilities 

Solve  for  the  partial  currents  using  the  face-to-face  transport  probabilities 
Compute  the  within-group  scattering  source  using  the  partial  current  solution 
and  the  angular  flux  distributions 
Refine  the  angular  flux  distributions 
until  the  angular  flux  distributions  are  converged 

Use  the  angular  flux  distributions  and  the  partial  currents  to  calculate  the  angular 
and  scalar  fluxes. 

Algorithm  1.1:  Overview  of  the  distribution  iteration 
algorithm. 


1 . 2  Mo  tivation 

A  reliable  and  robust  time-dependent  neutron  transport  code  that  can  rapidly 
converge  to  a  solution  is  useful  for  a  variety  of  applications,  one  example  being  the 
determination  of  activation  products  in  a  fast-flux  reactor.  All  current  discrete  ordi¬ 
nates  implementations  of  time-dependent  transport  use  the  source  iteration  algorithm 
and  exhibit  poor  performance  for  high  scattering  problems  (Lewis  &  Miller  [22]  and 
Reed  [29]). 

Prior  to  my  research,  the  distribution  iteration  algorithm  had  been  implemented 
in  slab  and  XY  geometry  for  monoenergetic  problems  and  isotropic  emission  sources. 


For  distribution  iteration  to  be  useful  in  a  production  code,  it  needed  to  be  extended 
to  support 

•  XYZ  geometry 

•  multigroup  transport 

•  time  dependence 

To  have  confidence  that  the  distribution  iteration  algorithm  is  a  suitable  replacement 
for  accelerated  source  iteration,  I  needed  to: 

•  Implement  a  formal  test  plan  that  demonstrates  that  the  distribution  iteration 
algorithm  correctly  and  reliably  solves  the  discrete  ordinates  equations. 

•  To  demonstrate  the  performance  advantage  that  distribution  iteration  has  over 
accelerated  source  iteration. 

1.3  Statement  of  the  Problem 

The  problem  was  to  develop  a  robust  time-dependent  algorithm  using  the  distri¬ 
bution  iteration  methodology  and  compare  the  performance  relative  to  an  established 
code,  e.g.  PARTISN.  The  extension  from  slab  to  XY  geometry  is  complicated  because: 

•  The  angular  quadratures  are  more  complicated. 

•  The  implementation  of  spatial  quadratures,  particularly  a  characteristic  method, 
is  more  complicated. 

•  Mesh  sweeps  involve  more  than  one  pair  of  faces. 

•  The  sparse  matrix  structure  used  to  solve  the  partial  current  problem  has  a 
more  complex  structure. 

The  extension  from  XY  to  XYZ  geometry  is  conceptually  simpler  than  the  extension 
from  slab  to  XY  geometry  because  the  angular  quadratures  are  similar,  the  spatial 
quadratures  are  derived  in  a  similar  fashion,  and  the  complexity  of  the  mesh  sweeps 
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and  the  sparse  matrix  are  comparable.  However,  a  “just  one  more  axis”  approach 
based  on  the  work  done  by  Prins  [28]  would  make  a  complicated  implementation  into 
an  unmanageable  software  engineering  problem  and,  hence,  difficult  to  verify. 

The  DI  algorithm  as  developed  by  Wager  and  Prins  sometimes  failed  to  con¬ 
verge  when  only  one  or  two  inner  iterations  were  used  to  refine  the  angular  current 
distribution  per  call  to  the  partial  current  solver.  Prins  was  able  to  achieve  conver¬ 
gence  for  his  test  problems  by  performing  up  to  ten  refinement  iterations  per  call 
to  the  partial  current  solver.  Though  this  method  did  improve  convergence,  it  was 
not  clear  that  it  would  work  in  all  cases  and  it  was  an  inefficient  approach.  Prins 
attempted  several  other  approaches  to  improve  the  refinement  of  the  angular  current 
distribution,  however,  none  of  them  were  completely  successful. 

1.4  Goals  of  the  Research 

The  work  done  by  Wager  [32]  and  Prins  [28]  demonstrated  the  robust  perfor¬ 
mance  of  DI  for  time-independent  problems  in  slab  and  XY  geometries.  The  conven¬ 
tional  algorithms  for  solving  multiple  energy  groups  and  time  dependent  problems 
work  by  modifying  the  emission  density,  q ,  and  the  total  cross  section,  a.  The  use  of 
DI  with  these  algorithms  has  not  been  previously  demonstrated  and  it  is  uncertain 
whether  the  changes  to  the  emission  density  and  the  total  cross  section  will  adversely 
impact  the  performance  of  the  DI  algorithm. 

The  goals  of  this  research  were: 

•  Solving  the  time-dependent  BTE  using  DI. 

•  Extending  the  time-dependent  DI  methodology  to  multiple  energy  groups. 

•  Extending  DI  to  support  XYZ  geometry  in  a  straightforward  and  methodical 
manner. 

•  Improving  the  performance  of  the  angular  current  distribution  refinement  algo¬ 
rithm. 
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•  Demonstrating  that  DI  is  a  viable  substitute  for  source  iteration. 

The  ultimate  goal  of  my  research  was  to  build  upon  the  previous  work  and  demonstrate 
that  DI  is  a  viable  alternative  to  source  iteration,  either  unaccelerated  or  accelerated, 
for  a  broad  range  of  transport  problems. 

1 . 5  Scope 

To  support  the  goals  of  my  research,  a  multigroup,  time-dependent  transport 
code  utilizing  DI  was  created.  The  code  supports  slab,  XY  and  XYZ  geometries 
using  a  variety  of  spatial  and  angular  quadratures.  Time-dependent  external  sources 
are  supported  and  there  is  limited  support  for  time-dependent  material  properties. 
The  formal  verification  effort  was  focused  on  slab  geometry  and  XY  geometry.  The 
verification  of  the  XYZ  geometry  implementation  was  designed  to  take  advantage  of 
common  design  basis  between  the  XY  and  XYZ  geometry  implementations. 

The  verification  testing  was  designed  with  two  separate  requirements.  The  first 
requirement  was  to  exercise  specific  aspects  of  the  DI  algorithm  in  order  to  check  for 
error  in  the  corresponding  code  segments.  The  second  requirement  was  to  test  the 
parameter  space,  e.g.  scattering  ratio,  in  order  to  demonstrate  that  DI  correctly  and 
robustly  solves  the  discrete  ordinates  equations. 

The  data  collected  during  the  verification  testing  was  then  used  to  evaluate  the 
performance  of  both  the  distribution  iteration  and  source  iteration  algorithms. 

1.6  Assumptions  and  Limitations 

Isotropic  cross-sections  are  used  in  order  to  focus  the  research  on  the  implemen¬ 
tation  of  the  DI  method  to  time-dependent  transport.  While  implementing  higher- 
order  terms  to  support  anisotropic  scatter  would  not  be  difficult,  it  does  add  com¬ 
plexity  to  the  software  and  complicate  verification  testing. 

The  time-dependent  material  properties,  such  as  density  and  composition,  were 
held  constant  during  testing.  The  research  did  not  address  the  implicit  relationship 
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between  neutron  flux  and  material  properties.  The  software,  however,  was  designed 
to  allow  for  changes  in  material  properties. 

The  multigroup  energy  support  in  my  research  was  limited  to  non-fissioning 
materials.  While  adding  support  for  fissioning  problems  is  not  complicated  algorith¬ 
mically,  the  effort  required  was  better  utilized  towards  the  other  research  goals.  A 
demonstration  of  multigroup  support  based  on  up  and  downscatter  is  more  than  suf¬ 
ficient  to  show  the  applicability  of  distribution  iteration  as  a  replacement  for  source 
iteration. 

The  performance  metric  was  limited  to  iteration  count  rather  than  run  time 
because  iteration  count  provides  an  indication  of  the  efficiency  of  the  code.  Run 
time  is  a  subjective  assessment  of  a  code  because  design  requirements  may  result  in 
different  performance  trade-offs. 

1 . 7  Approach 

The  first  step  was  to  determine  a  mathematical  formulation  of  the  time-dependent 
BTE  that  was  amenable  to  the  distribution  iteration  approach  and  maintains  a  close 
connection  to  the  physics.  Next,  a  one-dimensional,  time-dependent  distribution  iter¬ 
ation  was  implemented  for  monoenergetic  problems  and  the  performance  of  that  code 
was  compared  relative  to  PARTISN.  Finally,  the  development  of  metrics  for  compar¬ 
ing  the  performance  relative  to  a  benchmark  was  undertaken.  Some  of  the  points  of 
investigation  were  to: 

•  Evaluate  the  application  of  more  robust  spatial  quadratures,  such  as  first  mo¬ 
ment  and  characteristic  methods. 

•  Evaluate  time  stepping  methods. 

•  Evaluate  methods  to  improve  the  propagation  of  global  neutron  partial  currents. 

The  first  version  of  the  time-dependent  DI  code  (DI-1)  highlighted  the  require¬ 
ment  for  a  more  effective  algorithm  for  propagating  global  neutron  partial  currents  in 
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deep-penetration  problems.  Addressing  that  issue  was  undertaken  during  the  devel¬ 
opment  of  the  XY  geometry  version  of  the  code  (DI-2). 

The  implementation  of  an  improved  time-stepping  method  was  the  next  phase 
of  the  research.  The  Dl-1  code  developed  to  support  my  research  prospectus  was 
implemented  using  an  implicit  backward  difference  method  to  numerically  approxi¬ 
mate  the  time  derivative.  While  this  approach  was  straightforward,  truncation  error 
affected  the  performance.  The  most  obvious  approach  was  to  implement  the  mid¬ 
point  average  method  in  lieu  of  the  backward  difference  method  utilized  by  Dl-1.  The 
second  version,  Dl-2,  implemented  both  time  stepping  methods. 

The  implementation  of  XYZ  geometry  was  accomplished  in  the  third  version 
(Dl-3)  of  my  transport  code.  The  final  phase  was  the  extension  of  the  code  to  support 
multiple  energy  groups.  This  was  the  last  phase  because  it  is  dependent  on  having  a 
robust  monoenergetic  implementation. 
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II.  Theory 


Numerical  methods  for  solving  time-dependent  (evolutionary)  partial  differen¬ 
tial  equations  can  be  characterized  by  how  the  time  variable  is  treated  and  how  the 
extrapolation  to  a  future  state  is  performed.  If  all  the  derivatives-both  spatial  and 
time-are  represented  with  finite  differences,  then  the  method  is  fully  discretized  (Iser- 
les  [17]).  For  example,  given 


du  d  2u 
dt  dx 2 ’ 


(2.1) 


then  a  fully  discretized  scheme  using  the  forward  Euler  method  can  be  expressed  as 


w 


j+i 


=  uJ,  + 


At 

(Ax) 


(«;_!  -  2“f + “J+i) . 


(2.2) 


where  At  /(Ax)2  is  the  Conrant  number,  and  u  is  a  vector  of  spatial  values  indexed  by 
l  at  times  j  and  j  +  1.  An  alternative  method  is  the  semi-discretization  method.  In 
this  formulation  a  partial  differential  equation  is  transformed  into  a  coupled  system 
of  ordinary  differential  equations.  My  research  is  based  on  fully  discretized  methods. 

The  means  of  extrapolating  to  the  next  time  step  can  be  categorized  as  either 
explicit  or  implicit.  An  explicit  method,  such  as  Euler’s  method,  only  utilizes  values 
from  the  current  time  step  to  compute  values  for  the  next  time  step.  For  example, 
given 

Ait 

^  =  /(M)  (2.3) 

the  time  derivative  is  approximated  by  using  a  truncated  Taylor  series  expansion  as 


d  u  uj+1  —  uj 

dt  p  At 


(2.4) 


The  differential  equation  can  be  solved  for  each  time  step  by 


vP+1  =  v?  +  A tf(u\  f). 


(2.5) 
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Because  the  error  term  in  the  truncated  Taylor  series  expansion  is  0(Af2),  the  global 
truncation  error  is  first  order.  While  the  forward  Euler  method  is  convergent  (Iserles 
[17]),  it  is  only  conditionally  stable.  An  alternative  explicit  method  is  the  explicit 
midpoint  method  defined  as 

uj+2  =  uj  +  2  A tf(uj+1,tj+1).  (2.6) 


The  explicit  midpoint  method  is  classified  as  a  multistep  method  and  is  second  order. 

An  implicit  method,  such  as  the  trapezoidal  rule,  utilizes  values  from  multiple 
time  steps,  including  extrapolated  values,  to  compute  values  for  the  next  time  step. 
For  example,  for  the  backward  Euler  method,  the  time  derivative  is  approximated  as 


d  u  uj+1  —  uj 

d t  P+1  At 


(2.7) 


Thus,  the  differential  equation  is  solved  by 


uj+1  =  uj  +  A  tf(uj+1,tj+1). 


(2.8) 


The  backward  Euler  method  has  the  same  local  truncation  error  as  the  forward  Euler 
method  and  is  also  convergent.  Unlike  the  forward  Euler  method,  the  backward  Euler 
method  is  unconditionally  stable,  which  is  a  significant  advantage  when  dealing  with 
stiff  problems  (Iserles  [17]). 

The  primary  disadvantage  to  using  the  backward  Euler  method  is  the  first  order 
global  truncation  error.  The  trapezoidal  rule, 

uj+i  =  u3  +  ^  ( f(uJ+1 , tJ+1)  -  ,  (2.9) 

has  second  order  global  truncation  error  (Iserles  [17]),  is  convergent  and  is  uncondi¬ 
tionally  stable.  The  disadvantage  to  using  the  trapezoidal  rule  to  solve  the  BTE  is  the 
additional  memory  requirement  for  storing  /(i/7+l , U+I)  and  /(id, P).  The  implicit 
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midpoint  method  evaluates  /  at  the  midpoint  between  the  steps  j  and  j  + 1  and  yields 


+  (2.10) 

and  also  has  second  order  local  truncation  error. 

An  implicit  method  can  be  solved  by  making  it  into  a  semi-implicit  method, 
which  means  the  method  has  been  solved  by  linearization.  The  strategy  behind 
a  semi-implicit  method  is  to  split  the  problem  into  two  portions,  one  that  will  be 
solved  implicitly  and  another  that  will  be  solved  explicitly.  An  example  of  a  semi- 
implicit  method  is  PARTISN’s  approach  to  solving  the  time-dependent  BTE,  which 
is  discussed  in  detail  in  section  2.4.4  and  Appendix  A.  A  simple  illustration  of  a 
semi-implicit  method  can  be  shown  by  solving  the  backward  Euler  method  for  /, 

7  / 1 V 1  —  ti.j 

Hu1+1,t>+1)=  —  ■  (2.11) 

Substituting  for  /  into  the  forward  Euler  method  yields 

uj+ 2  =  2 uj+1  -  uj.  (2.12) 

Performing  a  Taylor  series  expansion  of  tP+2  and  tP  about  tP+1  yields 

uj+ 2  =  uj+1  +  A  tf(uj+1,tj+1)  +  tj+1)  +  0((At)3)  (2.13) 

and 

uj  =  uj+1  -  A tf(uj+1,  tj+1)  +  tj+1)  -  0((At)3).  (2.14) 

Taking  the  difference  between  uJ+2  and  uJ  yields 

uj+ 2  -  uj  =  2A tf(uj+1,ti+1)  +  0((At)3),  (2.15) 

which  is  the  explicit  midpoint  method  and  is  second  order. 


16 


Before  delving  into  the  implementation  of  the  time-dependent  BTE  in  DI,  it  is 
useful  to  cover  some  background  involving  the  time- independent  BTE. 

2.1  Operator  Notation 

For  clarity,  the  BTE  can  be  represented  in  operator  form  as 

d 

V1  —  T  =  -5T  +  Qe  (2.16) 

where  V  is  a  diagonal  matrix  of  neutron  group  velocities,  T  is  a  vector  containing 
angular  fluxes  for  each  group,  and  Qe  is  a  vector  of  embedded  sources.  The  transport 
operator  B  is 


BT  =  (L-S)V 

(2.17) 

(L'V)g  =  Cl  •  +  cr t(r,  t)^g 

(2.18) 

(S*)g  =  f  dCl'  Kg^g(r,  Cl'  -  0)vlv(f,  Cl'), 

9' 

(2.19) 

where  K  represents  the  scattering  kernel  and  at  is  the  total  cross  section.  The  operator 
L  represents  losses  and  the  operator  S  represents  scattering  sources. 

2. 2  Spatial  Quadratures 

Many  spatial  quadratures  have  been  developed  for  use  with  the  discrete  ordi¬ 
nates  formulation  of  the  BTE.  The  spatial  quadratures  can  be  broadly  grouped  as 
zeroth  moment  and  first  moment  methods.  The  zeroth  moment  methods  only  use  cell 
average  values,  e.g.,  cell  average  angular  fluxes,  while  the  first  moment  methods  use 
both  the  cell  average  values  and  the  first  spatial  Legendre  moments. 

The  most  common  zeroth  moment  methods  are  the  Diamond  Difference  (DD), 
Step  Characteristic  (SC)  (Lathrop  [20]),  and  Step  approximations.  The  diamond 
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difference  spatial  quadrature  uses  the  auxiliary  relationship  (in  slab  geometry) 

^l  +  ^r  =  2^a,  (2.20) 

where  and  are  the  angular  fluxes  on  the  left  and  right  boundaries,  respectively. 
The  advantage  to  using  diamond  difference  is  the  simplicity  of  its  implementation 
and  second-order  accuracy.  The  primary  disadvantage  is  the  fact  that  it  will  gener¬ 
ate  physically  meaningless  negative  fluxes  in  some  circumstances.  The  step  spatial 
quadrature  uses  the  auxiliary  relationships 

^R  =  ^A  p  >  0  (2.21) 

p  <  0.  (2.22) 

The  step  spatial  quadrature  is  simpler  than  diamond  difference,  however,  it  only  has 

first-order  convergence,  which  makes  it  impractical  to  use  in  a  production  code.  It  is, 
however,  a  “positive”  method  in  the  sense  that  it  will  always  generate  positive  fluxes. 
The  step  characteristic  spatial  quadrature  is  based  on  the  assumption  that  particles 
travel  in  straight  lines  between  collisions.  These  straight  lines  are  the  characteristics 
of  the  BTE  and  are  represented  by 

=  q,  (2.23) 

where  s  is  the  path  traveled  by  the  particle.  Without  loss  of  generality,  consider  the 
p  >  0  case.  Integrating  the  characteristic  equation  yields 

ip(s)  =  ijjLe~atS  +  f  d s'  q(s,)e~(Tt('s~s'\  (2.24) 

Jo 

The  step  characteristic  spatial  quadrature  has  the  simplifying  approximation  that 
q(s)  =  qA,  where  qA  is  the  cell  average  source.  The  advantage  to  using  SC  is  that  it 
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is  a  positive  method  and  has  second-order  convergence.  The  drawback  to  SC  is  that 
it  is  difficult  to  generalize  to  non-rectangular  coordinate  systems. 

The  most  common  first  moment  methods  are  Linear  Discontinuous  (LD)  (Hill 
[15])  and  Linear  Characteristic  (LC)  (Larsen  &  Alcouffe  [18]).  The  first  moment 
methods  have  the  drawback  of  increased  memory  requirements,  however,  they  do 
exhibit  a  better  rate  of  convergence  than  the  zeroth  moment  methods.  The  linear 
discontinuous  spatial  quadrature  assumes  that  the  flux  is  linear  within  a  cell  but 
discontinuous  with  the  flux  in  adjacent  upstream  cells.  This  method  can  produce 
negative  fluxes  with  optically  thick  cells.  The  linear  characteristic  method  is  similar 
to  the  step  characteristic  method  with  a  first-moment  expansion  for  the  source 

q(x)  =  qAP0(x )  +  qxPi(x),  (2.25) 

where  qx  is  the  first  moment  and  Pn  are  the  Legendre  polynomials  shifted  and  scaled 
in  x.  Though  LD  is  not  a  positive  method,  it  produces  fewer  negative  fluxes  than 
DD. 

There  are  new  spatial  quadratures  that  have  been  developed  which  have  not 
gained  widespread  acceptance.  An  important  reason  for  the  lack  of  acceptance  is  due 
to  the  complexity  of  developing  an  accelerator  for  use  with  source  iteration.  Some 
examples  of  the  newer  spatial  quadratures  are  Exponential  Discontinuous  (Wareing 
[34]),  Nonlinear  characteristic  (Walters  &  Wareing  [33]),  Exponential  Characteristic 
(Mathews,  Sjoden  &  Minor  [25])  and  Nonlinear  Corner  Balance  (Castrianni  &  Adams 
[10]) 

The  negative  fluxes  can  have  a  deleterious  effect  on  the  performance  of  a  trans¬ 
port  code  (Lathrop  [20]).  In  multidimensional  problems,  the  negative  fluxes  can  slow 
the  rate  of  convergence  or  render  the  acceleration  method  used  with  source  iteration 
ineffective.  The  propagation  of  negative  fluxes  can  result  in  a  lack  of  stability  for  time 
dependent  problems.  Traditionally,  a  negative  flux  fix-up  is  typically  used  with  non- 
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positive  methods  like  DD.  The  drawback  to  using  a  negative  flux  fix-up  is  a  slower 
rate  of  convergence. 

2.3  Distribution  Iteration 

Distribution  iteration  solves  the  BTE  by  solving  two  different  aspects:  the  par¬ 
tial  current  problem  and  the  distribution  of  the  angular  partial  currents  on  the  cell 
faces.  The  partial  currents  problem  is  a  simplification  of  the  transport  problem  into  a 
linear  algebra  problem  that  can  be  solved  directly.  The  second  aspect,  solving  for  the 
distribution  of  the  angular  partial  currents,  serves  two  purposes-providing  the  means 
for  improving  the  partial  currents  problem  and,  second,  providing  the  solution  to  the 
angular  and  scalar  fluxes  in  each  cell. 

Before  proceeding,  a  brief  explanation  on  notation  will  help  guide  the  reader 
through  the  derivations  that  follow.  The  use  of  the  vector  decoration,  e.g.,  J,  denotes  a 
vector  constructed  from  multiple  faces.  Symbols  which  are  inherently  vectors  because 
of  the  use  of  an  angular  quadrature,  such  as  the  cell  average  angular  flux  (-0A),  are 
not  decorated.  The  rationale  for  restricting  the  use  of  the  vector  decoration  was  to 
make  it  clear  to  the  reader  which  variables  are  constructed  from  multiple  faces.  The 
use  of  the  “blackboard  bold”  style,  e.g.  M,  denotes  a  matrix.  I  have  also  adopted 
the  convention  of  using  the  “e”  superscript  to  denote  all  external  contributions  to  the 
within-group  BTE. 

The  DI  algorithm  models  the  BTE  as  a  collection  of  cells  coupled  together  by 
partial  current  flows  of  neutrons  between  the  cells.  An  implicit  representation  of  the 
global  partial  currents  problem  can  be  expressed  as  (Mathews  [24]) 

J  =  Je  +  MJ,  (2.26) 

where  J  is  the  vector  form  of  all  the  partial  currents  across  cell  faces;  Je  is  the  emission 
partial  current  due  to  scatter  from  other  groups,  fission,  and  fixed  emitters;  and  M  is 
a  matrix  of  the  mean  face-to-face  transport  probabilities. 
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The  advantage  of  (2.26)  is  that  J  is  non-negative  because  each  element  of  M  is 
non-negative,  as  is  each  element  of  Je.  The  implicit  form  does  have  the  disadvantage 
that  an  iterative  solution  will  exhibit  slow  convergence  for  optically  thick  cells  as 
c  — »  1  (Mathews  [24]).  Solving  for  J  explicitly  yields 

J  —  (I—  M)_1  Je,  (2.27) 

where  I  is  the  identity  matrix.  If  we  apply  synthetic  division  to  the  matrix  inversion, 
we  have 

(I  —  M)-1  =  1  +  M  +  M2  +  ...,  (2.28) 

thus  the  inverted  matrix,  which  we  identify  as  G,  represents  all  possible  transport 
paths  that  a  neutron  might  take  over  the  entire  spatial  domain. 

A  method  for  determining  the  elements  of  the  matrix  M  is  required.  The  first 
step  is  to  formulate  the  equations  for  within-cell  transport1  as 

+  Kv’sS^a  (2.29) 

and 

Jout  =  KOIJm  +  KOE?o  +  KOSE^A;  (2.30) 

where  ,iJja  is  the  cell  average  angular  flux,  jm  is  the  inward  flowing  currents  organized 
on  a  per-cell  basis  and  jout  is  the  outward  flowing  currents  organized  on  a  per-ccll 
basis.  The  physical  significance  of  each  matrix  is 

•  K01  accounts  for  the  contribution  from  the  uncollided  neutrons  to  the  cell  out¬ 
flow, 

•  Koe  accounts  for  the  contribution  from  the  cell  fixed  emission  source  to  the  cell 
outflow, 

T  have  adopted  an  angular  partial  current  notation  vice  an  angular  flux  notation. 
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•  IKos  accounts  for  the  contribution  from  the  cell  scattering  source  to  the  cell 
outflow, 

•  KA1  accounts  for  the  contribution  from  the  nncollided  neutrons  to  the  cell  aver¬ 
age  angular  flux, 

•  KAe  accounts  for  the  contribution  from  the  cell  fixed  emission  source  to  the  cell 
average  angular  flux, 

•  KAs  accounts  for  the  contribution  from  the  cell  scattering  source  to  the  cell 
average  angular  flux,  and 

•  E  accounts  for  the  scatters  within  the  cell. 

The  IK  matrices  are  dependent  upon  the  spatial  quadrature  scheme,  angular  quadra¬ 
ture,  and  material  properties. 

Solving  (2.29)  for  -0A  yields 

(I  _  k^sE)  iJja  =  K^ljm  +  K^V,  (2.31) 

where  I  is  the  identity  matrix.  Let 

L=  (I-.IK^E)-1.  (2.32) 

Applying  synthetic  division  to  the  matrix  inversion  yields 

L  =  I  +  KV’SE  +  (Kv'sE)2  +  ...  (2.33) 

Thus,  the  matrix  L  accounts  for  all  combinations  of  transport  within  a  ccll-uncollided, 
single  collision,  two  collisions,  etc.-and  the  cumulative  contribution  to  the  cell  average 
angular  flux.  Substituting  in  L  yields 

iJja  =  +  LK^V.  (2.34) 
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Substituting  (2.34)  into  (2.30)  yields 


Jout  =  (KOI  +  KOSSLKyi^)  Jm  +  ^KOE  +  KOSSLKyE^  ge.  (2.35) 

The  coefficient  matrix  and  external  angular  partial  current  vector  are  defined  as 

moi  =  K01  +  (2.36) 

and 

f  =  (K0E  +  KosELK^e)  qe.  (2.37) 

The  partial  currents  vector  Je  is  computed  from  the  angular  partial  currents  by  us¬ 
ing  the  angular  quadrature  weights.  The  angular  quadrature  weight  matrix,  W,  is 
constructed  such  that 

jout  =  wjout_  (2.38) 

The  equation  for  the  partial  current  for  each  cell  is  then 

J  =  WmoiJin  +  Wf .  (2.39) 


The  inward  flux  distribution,  £,  is  defined  as 


C faceiXi 


o  outward 

^  L  *  nface 

fr'e r/acedr/  /d'-n^tcyard<0 

r>/  4^  outward 

^  ‘  nface 

ip(r',  Cl') 

(2.40) 


where  is  the  ontward  normal  on  a  cell  face  face.  The  flnx  distribution  ( 

apportions  the  partial  current  on  a  cell  face  to  the  individual  ordinates  that  constitute 
the  inward  flow  into  a  cell.  The  inward  flnx  distribution  can  be  represented  in  matrix 
form  as 

7in 

ryi  _  J  cell,  face 

^ cell, face  ~T 

cell, face 


(2.41) 
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where  Z ceujace  is  a  subset  of  Z  that,  for  a  given  cell  and  face,  j™u.face  is  the  portion  of 
jm  for  a  cell  that  corresponds  to  face  and  Jceii,face  is  the  inward  partial  current  for  a 
cell  that  corresponds  to  face,  which  is  a  scalar  value.  Using  Z  to  eliminate  jm  yields 

J  =  WmoiZ  J  +  Wj°ut.  (2.42) 

Thus,  we  can  determine  M  from 

M  =  WmoiZ  (2.43) 

and  Jc  from 

Je  =  W f.  (2.44) 

The  basic  algorithm  for  DI  that  1  implemented  in  my  research  is  summarized 
in  Algorithm  2.1.  There  are  two  key  points  in  the  algorithm  that  differ  from  the 
implementations  by  Wager  and  Prins  that  I  will  briefly  highlight.  The  most  significant 
point  of  departure  is  in  the  sweeping  methodology.  The  cell  outward  angular  partial 
currents  can  be  determined  from  the  equation 

Jout  _  KOIJm  +  KOSELKyiJin  +  Je.  (2.45) 

The  cell  outflow  due  to  within-cell  sources  is  defined  as 

r  =  A'r+f,  (2.46) 

where 

A'  =  KosSLK^’1.  (2.47) 

This  change  to  the  algorithm  solves  the  stability  problem  that  was  inherent  in  the 
previous  algorithm.  The  algorithm  used  by  Wager  and  Prins  would  update  the  scat¬ 
tering  source  while  sweeping  while  this  algorithm  computes  the  scattering  source  for 
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all  cells  before  sweeping.  The  Wager  and  Prins  algorithm  uses  partially  updated  cell 
inflows  when  calculating  the  scattering  source,  which  can  cause  unstable  behavior, 
particularly  in  XY  and  XYZ  geometry.  My  algorithm  also  has  improved  parallel  ef¬ 
ficiency  because  each  cell  can  be  calculated  independently  of  the  others.  Also,  the 
contribution  due  to  sources  in  each  cell  can  be  computed  for  all  ordinates  in  a  cell 
using  a  sparse  BLAS  library. 

The  other  point  that  bears  further  explanation  is  the  phrase  “Until  Z,  stabilizes.” 
Typically,  only  one  iteration  is  needed  to  refine  the  inward  current  distribution.  Dur¬ 
ing  the  course  of  my  research,  there  were  some  configurations  that  required  more  than 
one  iteration  in  order  for  DI  to  converge,  but  two  always  sufficed  (in  my  testing).  Al¬ 
ternatively,  one  can  iterate  until  the  inward  current  distribution  converges.  In  fact, 
it  is  possible  to  solve  the  transport  equation  this  way  without  explicitly  solving  for 
the  partial  currents,  but  the  “no  partial  current  solver”  is  a  poor  algorithm:  The 
overall  number  of  iterations  increases  because  the  refinement  of  the  inward  current 
distribution  does  not  propagate  the  cell-to-cell  transport  efficiently.  There  is,  how¬ 
ever,  an  optimization  between  the  number  of  iterations  required  to  refine  the  inward 
current  distribution  in  order  to  converge  and  the  relative  CPU  cost  between  solving 
the  partial  currents  problem  and  refining  the  inward  current  distribution. 

Because  the  computation  of  the  angular  flux  distribution  within  the  cell  can 
be  performed  independently  of  any  other  cell,  parallelization  can  be  accomplished 
with  modest  effort.  Furthermore,  parallel  implementations  of  linear  solvers  are  well 
established,  which  benefits  the  solution  of  the  global  partial  currents  problem. 

2.3.1  Partial  Currents  in  Slab  Geometry.  Before  discussing  the  partial 
currents  in  XYZ  geometry,  an  example  from  slab  geometry  illustrates  the  basic  termi¬ 
nology.  Figure  2.1  is  a  simple  example  illustrating  the  coupling  between  two  slabs.  Let 
Jf  and  J~  represent  the  partial  currents  in  the  positive  /x  and  negative  [i  directions 
for  face  i,  respectively. 
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Compute  transport  matrices  (e.g.  IK,  E,  and  W) 

Initialize  inward  current  distribution,  Z  (e.g.  isotropic) 

Compute  within  group  source  (je  and  Je) 
repeat 

Compute  M 

Solve  (I  -  M)  J  =  Je 

repeat 

Compute  the  source  outflow  in  each  cell  as  jos  =  je  +  A'Z  J 
for  each  octant  do 

Initialize  inward  angular  current,  jm  from  boundary 
for  each  cell  in  flow  sequence  do 

Compute  outward  cell  angular  partial  currents  jout  =  IKOIjm  +  jos 
Propagate  outward  cell  angular  partial  currents  to  neighboring  cells  as 
new  inward  cell  angular  partial  currents 

end  for 
end  for 

Update  Z  using  updated  inward  cell  angular  partial  currents 
until  Z  stabilizes 
until  Z  converges 

Algorithm  2.1:  Distribution  Iteration  Algorithm 


Figure  2.1:  Coupled  Cells  in  Slab  Geometry 
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The  partial  currents  shown  in  Figure  2.1  can  be  expressed  as 


4- 
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i 

(2.48) 
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using  the  notation  to  represent  the  mean  probability  that  a  neutron  will 

travel  from  face  f  rom  to  face  to  and  J^e  is  the  fixed  source  current  on  a  face  in 
either  a  positive  or  negative  /i  direction.  When  constructing  the  current  vector  J 
from  the  individual  4  ■  I  adopted  the  following  convention: 


J  = 


4 

4 

4 

4 


(2.54) 


The  individual  m’s  are  used  to  construct  a  linear  system 
vector,  specifically  (I  —  M)  J  =  Je,  where 


0  MfA 
0 


0 

0 


0 

M+ 


that  solves  for  the  J  current 


0 

0 

0 

0 


(2.55) 
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2.3.2  Partial  Currents  in  XYZ  Geometry.  Prins  [28]  developed  an  XY 
geometry  implementation  of  DI.  Conceptually,  the  partial  currents  problem  in  XY 
geometry  can  be  considered  the  same  as  four  slab  geometry  partial  currents.  The 
four  arrangements  are  the  Left/Right  coupling,  the  Top/Bottom  coupling,  and  the 
two  diagonal  couplings,  Left/Bottom  to  Right/Top  and  Right/Bottom  to  Left/Top. 
Extending  into  XYZ  geometry  involves  extending  those  couplings  into  the  Z  direction. 

By  deconstructing  the  partial  currents  and  treating  all  the  cell  faces  identically,  a 
simple  arrangement  of  cell  faces  into  one  vector  is  possible.  One  of  my  implementation 
goals  was  to  maximize  the  common  code  elements  between  the  different  geometries. 
That  goal  leads  to  an  extensible  arrangement  of  the  partial  currents  vector  where  all 
the  faces  normal  to  the  X  axis  (X  faces)  are  grouped  together,  all  the  faces  normal 
to  the  Y  axis  (Y  faces)  are  grouped  together,  and  all  the  faces  normal  to  the  Z  axis 
(Z  faces)  are  grouped  together.  In  a  slab  geometry  problem  of  Nx  cells,  there  are 
Nx  +  1  faces.  In  XY  geometry  of  Nx  by  Ny  cells,  there  will  be  (Nx  +  l)iVy  X  faces 
and  (Ny  +  1  )NX  Y  faces.  The  partial  current  vector  would  then  consist  of  a  vector 
containing  a  block  of  X  faces  and  block  of  Y  faces.  In  XYZ  geometry  of  Nx  by  Ny  by 
Nz  cells,  there  will  be  (Nx  +  l)NyNz  X  faces,  (Ny  +  1)NXNZ  Y  faces  and  (Nz  +  l)iVxIVy 
Z  faces.  The  partial  current  vector  would  be  extended  by  a  block  of  Z  faces. 

This  arrangement  of  partial  currents  leads  to  an  interesting  arrangement  of  the 
matrix  I  —  M.  The  terms  adjacent  to  the  main  diagonal  are  associated  with  the  axial 
Left/Right,  Bottom/Top,  and  Back/Front  transport.  Off  the  diagonal  are  clusters 
that  are  associated  with  the  diagonal  transport,  e.g.,  Left/Top.  Figure  2.2  shows  a 
graphical  representation  of  the  G  matrix.  The  columns  arc  labeled  with  the  associated 
inflow  dimension  and  the  rows  are  labeled  with  the  associated  outflow  dimension. 

2.3.3  Transport  Coefficients  in  XYZ  Geometry.  Appendix  B  shows  the 
derivation  of  the  transport  coefficients  for  the  step  and  diamond-difference  spatial 
quadratures  in  XYZ  geometry.  For  diamond-difference  to  have  positive  fluxes,  all  the 
transport  coefficients  must  be  positive.  All  the  coefficients,  except  for  the  three  axial 


Figure  2.2:  Graphical  Representation  of  the  Global  Par¬ 
tial  Currents  Matrix 


transmission  coefficients,  are  unconditionally  positive.  The  three  axial  transmission 
coefficients,  IK®*,  and  IK®1,  are  positive  if 
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respectively,  are  satisfied. 

However, 

it  is 

impossible  for  all  three  to  hold. 

In  fact,  at 

most  only  one  can  be  satisfied.  In  the  S„  angular  quadrature,  the  positive  coefficients 
will  be  the  ones  closest  to  the  axis.  For  example,  if  p  >  rj  +  £  for  an  ordinate,  its  IK®^ 
is  positive  and  its  K®y:  and  IK®1  are  negative. 

The  consequence  of  negative  coefficients  is  that  the  angular  flux,  in  the  absence 
of  sufficiently  strong  within-cell  sources,  can  be  biased  towards  axial  flows.  This  effect 
is  not  unique  to  DI,  it  will  also  occur  in  the  source  iteration  algorithm.  Thus,  in  XYZ 
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geometry,  the  use  of  diamond-difference  can  result  in  anisotropic  flux  distributions  as 
an  artifact  of  the  spatial  quadrature. 

In  the  DI  algorithm,  negative  fluxes  can  cause  problems  with  the  inward  angular 
flux  distribution.  If  some  of  the  elements  of  jm  are  negative  and  the  inward  partial 
current  remains  positive,  then  (  will  be  negative  where  jm  is  negative.  However,  if 
enough  elements  of  jm  are  negative  such  that  the  partial  current  is  negative,  then  ( 
will  have  the  opposite  polarity  with  respect  to  jm.  The  sign  of  (  also  affects  the  sign 
of  the  elements  of  M,  which  can  result  in  the  additional  production  of  negative  partial 
currents.  As  a  result,  the  elements  of  the  matrix  Z  can  oscillate  and  not  converge. 

The  conventional  method  for  implementing  a  negative  flux  fix-up  is  not  the  best 
approach  in  DI  because  it  requires  additional  storage  or  more  computations.  Instead, 
I  developed  a  straightforward  negative  flux  fix-up  scheme  that  takes  advantage  of  the 
normalization  of  (.  During  a  mesh  sweep,  if  a  negative  current  is  generated,  it  is  set 
to  zero.  When  the  matrix  Z  is  updated,  the  corrected  angular  partial  current  will  not 
contribute  to  the  partial  current,  thus  particle  conservation  is  assured. 

2.4  Time  Discretization 

The  most  established  time  dependent  BTE  transport  codes  are  the  1972  TIMEX 
code  (Reed  [30]),  the  1976  TIMEX  code  (Hill  [16])  and  PARTISN  (Alcouffe  [4]).  Some 
of  the  other  time-dependent  codes  that  have  been  developed  are  TDTORT  (Goluoglu 
[14]),  TDA  (Time-Dependent  ANISN)  (Dupree  [11]),  and  TRANZIT  (Lathrop  [21]). 

2-4-1  TIMEX  (1972)  Implementation.  The  TIMEX  (1972)  code  (Reed  [30]) 
utilizes  a  first  order  Taylor  series  expansion  of  the  time  derivative  to  discretize  the 
time  variable  analogous  to  forward  Euler.  TIMEX  (1972)  also  uses  the  multigroup 
approximation  to  discretize  energy  and  the  discrete  ordinates  representation  to  dis¬ 
cretize  in  angle.  The  BTE  solved  by  TIMEX  (1972)  in  slab  geometry  (Appendix  A) 
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in  operator  notation  is 


V_1-^'+1  +  W+1  =  SW  +  Qe.  (2.59) 

at 

Of  particular  interest  is  the  manner  in  which  the  transport  operator  B  from  (2.17) 
was  split.  The  loss  operator,  L,  operates  on  the  angular  flux  at  time  j  +  1  while  the 
scattering  source  operator,  S,  operates  on  the  angular  flux  at  time  j.  Reed  notes  [30] 
that  this  approach  is  inaccurate,  particularly  when  transient  events  occur,  e.g.,  a  delta 
distributed  source.  It  does  have  the  advantage  that  it  avoids  using  source  iteration 
to  solve  for  the  within-group  scattering  source.  The  TIMEX  1972  code  implements 
several  strategies  to  overcome  the  inherent  inaccuracy  of  this  approach  while  still 
avoiding  using  source  iteration  to  solve  for  the  within-group  scattering  source.  The 
details  of  those  strategies  are  not  relevant  to  my  research  because  DI  is  used  to  solve 
for  the  within-group  scattering  source. 

The  use  of  diamond  difference  as  a  spatial  quadrature  can  result  in  negative 
angular  fluxes.  As  shown  in  Appendix  A,  this  can  occur  when 

(i  +  ff*)>>0'  <2-60> 

hence  when  there  are  large  cells,  large  cross  sections,  or  short  time  steps.  TIMEX 
(1972)  has  a  zero  flux  fix-up  scheme  to  account  for  this  possibility. 

One  interesting  feature  used  by  TIMEX  (1972)  is  the  separation  of  the  first 
flight  flux  and  first  scatter  source  as  a  mechanism  to  handle  sources  delta  distributed 
in  space  and/or  time.  Delta  distributed  sources  are  problematic  for  any  transport 
code  because  short  time  steps  are  required  to  numerically  resolve  the  source.  This 
presents  two  problems:  first,  the  short  time  steps  can  cause  the  accumulation  of  error 
due  to  loss  of  precision;  and  second,  the  short  time  steps  can  result  in  negative  fluxes. 

Because  the  transport  operator  B  is  linear,  the  angular  flux  can  be  represented 
as  T  =  dT+Tc  where  T,,  is  the  uncollided  flux  and  Tc  is  the  collided  flux.  Substituting 
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into  (2.16)  yields 


v-1%-*»  +  l*u  =  q° 

at 

V-1^c  +  LVc  =  S(yu  +  yc 
at 


(2.61) 

(2.62) 


The  uncollided  flux  can  be  solved  analytically  and  its  first  collision  scattered  neutrons, 
SXVU,  are  then  used  as  the  fixed  source  when  solving  the  collided  flux  numerically. 

2-4-2  TIMEX  (1976)  Implementation.  The  TIMEX  (1976)  (Hill  [16])  code  is 
similar  to  the  TIMEX  (1972)  code,  the  key  differences  being  the  addition  of  delayed 
neutrons  and  the  introduction  of  a  synthetic  cross  section  and  a  synthetic  source. 
Consider  the  divided  difference  approximation 


d% b  ?/d+1  —  ?/d 

~dt  A t 


(2.63) 


where  ij)9+1  is  the  angular  flux  at  t/+i,  ^  is  the  angular  flux  at  tj,  and  At  is  the  time 
interval  t]+i  —  tj.  Substituting  (2.63)  into  the  BTE  yields 


VgAt 


(Vg+1  -ify+Cl-  V^+1  +  = 


g- 1 


G 


Y  +  Y  </  .A  +  <"■  (2-64) 


g'= i 


Collecting  the  ?/d+i  terms  and  moving  the  to  the  right  hand  side  yields 


Ci  •  v^0+1 


7+1  . 

(jf  T  — - — 
tg  vnA  t 


ti+1  = 


9- 1 


G 


+  Y  aYA  +  Y  +  irxiY  (2-65) 


g'= i 


VgAt  9 
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Define  the  synthetic  total  cross  section,  crt  ,  as 


aj+ 1  -  aj+1  H _ — 

°t q  ~  °t„  W  A  ,  ) 
9  9  Vn/Xt 


and  the  synthetic  source,  qg,  as 


■S+1  =  <4+1  + 


- ibj 

vnAV9 


(2.66) 


(2.67) 


so  that 


9- 1 


G 


•  W+1  +  a{+ liPi+1  = 


=  E- 


ri+! 


£+1 


+  E 


(j9 


+  ^+1, 


(2.68) 


g'= i  g'=g 

which  is  in  the  same  form  as  the  stationary  BTE.  The  TIMEX  (1976)  code  used 
a  time- independent  transport  code,  specifically  ONETRAN  (Hill  [15]),  to  solve  this 
stationary  BTE. 


The  authors  of  the  TIMEX  (1976)  (Hill  [16])  state  that  they  would  prefer  to 
solve  the  implicit  form 


V-1^-Tj+1  +  L¥+1  =  S^j+1  +  Qc 
at 


(2.69) 


The  implicit  form  was  not  used  because  of  the  amount  of  computer  time  that  would 
be  required.  For  example,  ONETRAN  required  5.5  minutes  on  a  CDC  7600  to  solve  a 
20  group,  134  interval  spatial  mesh,  S4  kes  calculation  (Hill  [15]).  The  transformation 
of  the  time-dependent  BTE  by  using  the  synthetic  cross  section  points  to  the  logi¬ 
cal  approach  of  using  an  established  time-independent  transport  code  to  solve  time 
dependent  problems. 


2-4-3  Implicit  Backward  Euler.  Before  proceeding  to  PARTISN,  it  is  useful 
to  examine  the  implicit  implementation  of  (2.69).  The  implicit  in  time  version  of  the 
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within-group  BTE  is 


Q  ■  W+1  +  a{^Lipi+i  =  a: 


j  ' J  -  1  - 


rs+1 

On—^a 


3+1 


3-1 

+  £ 

3'=1 


(7; 


7+1  +4+1 


G 

+  £ 

3,=3+l 


C7„ 


+#+1- 


(2.70) 


While  both  the  explicit  form,  which  is  used  by  TIMEX,  and  the  implicit  form  have  first 
order  local  truncation  error,  the  implicit  form  has  the  advantage  of  being  numerically 
stable  (Lewis  &  Miller  [22]). 

One  of  the  key  challenges  in  using  the  synthetic  cross  section  is  selecting  the 
time  step.  Consider  the  monoenergetic,  synthetic  BTE  in  slab  geometry  (suppressing 
the  indices  for  clarity) 


Ax 


ipL)  +  atipA  =  as(j)  +  q. 


(2.71) 


Let  e  be  the  optical  depth,  defined  as 


a  Ax 

e  = - 

/i 


(2.72) 


The  choice  of  spatial  quadrature  will  dictate  a  maximum  optical  thickness,  etoi,  that 
is  acceptable.  The  maximum  cell  size  is  then 


Ax- 


max  £ 


I  fin 


<?t 


"  £tol  • 


(2.73) 


Clearly  as  vAt  — >  0,  the  synthetic  cross  section  becomes  large  and  the  maximum  cell 
size  goes  to  zero.  Conversely,  as  vAt  — >  oo,  the  time  dependent  behavior  is  lost  and 
6t  -+•  o-t. 

The  Courant-Friedrichs-Lewy  condition  requires  that 


vAt 

Ax 


<  1. 


(2.74) 
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Making  the  substitution  for  Axmax  into  the  Courant-Friedrichs-Lewy  condition  yields 


vA  t 


iMmin  | 


<  1, 


(2.75) 


^tol 


which,  after  substituting  for  ert  to  get  crt,  yields  the  constraint 


0  <  At  <  -  ^l/imin|eto1  1 
o' t 


(2.76) 


Because  both  crt  and  v  are  positive,  then 


|  bin  in  |  Pol  1  b  0 


(2.77) 


must  be  true  in  order  to  have  positive  time  steps.  The  above  constraint  can  be  written 
as 

etoi  >  I — - — r,  (2.78) 

IP'min  | 

thus,  for  time  dependent  problems,  spatial  performance  decreases  as  the  angular 
quadrature  is  refined.  Picking  the  appropriate  time  step  requires  balancing  angu¬ 
lar  resolution  and  spatial  resolution  requirements.  This  constraint  also  makes  the  use 
of  non-positive  spatial  quadratures,  such  as  diamond-difference,  in  time-dependent 
transport  problematic.  Diamond-difference  in  slab  geometry,  for  example,  requires 
that 

Ax  <  2^^  (2.79) 

cr  t 

be  true  to  ensure  positive  fluxes.  Applying  the  requirement  for  Axmax  yields  the 
constraint 

etoi  <  2.  (2.80) 

Thus,  to  maintain  positivity  in  slab  geometry,  the  angular  quadrature  must  have 
/imin  >1/2,  which  in  the  Sn  angular  quadrature  is  only  true  with  S2.  If  a  negative 
flux  fix-up  scheme  is  implemented,  the  flux  will  be  biased  towards  the  ordinates  that 
meet  the  constraint. 


35 


While  it  is  possible  for  monoenergetic  problems  to  pick  a  time  step  that  results 
in  a  reasonable  optical  thickness,  e.g. 

At-—,  (2.81) 

vat 

the  speed  differential  between  thermal  and  fast  neutrons  makes  picking  a  uniform 
time  step  for  multigroup  problems  impossible.  Fast  neutrons  are  approximately  6000 
times  faster  than  thermal  neutrons,  which  can  result  in  synthetic  cross  sections  that 
are  10  to  1000  times  larger  between  energy  groups. 

Multigroup  problems  that  only  have  downscatter  can  mitigate  the  time  step 
issue  by  using  a  different  time  step  for  each  energy  group.  However,  if  there  is  both 
upscatter  and  downscatter,  this  approach  becomes  unworkable  because  temporal  res¬ 
olution  has  to  be  created  when  going  from  a  longer  time  step  energy  group  to  a  shorter 
time  step  energy  group.  Similarly,  if  different  spatial  grids  are  used  for  each  energy 
group,  spatial  resolution  has  to  be  created  when  going  from  a  coarse  grid  to  a  fine 
grid. 

Given  the  difficulty  of  picking  an  appropriate  time  step  that  works  equally  well 
for  all  energy  groups,  the  best  solution  would  be  to  use  a  spatial  quadrature  that  is 
non-negative  and  works  with  optically  thick  cells.  The  alternative  is  to  use  fixups 
in  the  transport  code  such  that  it  maintains  adequate  performance  with  optically 
thick  cells.  The  latter  approach  has  been  taken  in  production  codes  and,  thus,  D1 
has  to  demonstrate  equivalent  performance  with  optically  thick  cells  in  order  to  be  a 
substitute  for  source  iteration. 

2-4-4  PARTISN  Implementation.  PARTISN  implements  a  semi-implicit 
method  to  discretize  the  time  derivative  in  the  BTE.  Unlike  the  TIMEX  codes,  PAR¬ 
TISN  is  implicit  in  time  and  the  extrapolated  angular  fluxes  are  second  order  accurate. 
Appendix  A  shows  the  derivation  as  presented  by  Alcouffe  &  Baker. 
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Like  TIMEX  (1976),  PARTISN  utilizes  a  time-independent  code,  in  this  case 
DANTSYS,  to  solve  the  BTE  at  the  midpoint  between  time  steps.  If  the  extrapolated 
flux  is  negative,  a  set  to  zero  fix-up  strategy  is  utilized  in  conjunction  with  the  balance 
equation.  Because  PARTISN  uses  an  implicit  method  for  solving  the  midpoint  flux, 
it  is  numerically  stable.  The  set  to  zero  fix-up  can  result  in  less  than  second  order 
convergence. 

2-4-5  TDTORT  Implementation.  TDTORT  implements  a  flux  synthesis 
method  (Goluoglu  [14]),  which  assumes  that  time  dependence  is  partially  separable, 
i.e., 

d>(r,  ft,  E,  t )  =  T(f,  ft,  E,  t)T{t).  (2.82) 

The  function  $(r,  Cl,E,  t)  is  assumed  to  be  composed  of  a  shape  function,  T,  that  is 
slowly  varying  in  time  and  an  amplitude  function,  T,  that  is  typically  rapidly  varying 
in  time.  After  substituting  (2.82)  into  the  time-dependent  BTE  and  dividing  by  T(t), 
we  have 


ft  •  V  +  crt(r,  E ,  t)  T  (r,  ft,  E,  t )  = 

/oo  r 

d E'  /  dft'  crs (r ,  E'  ->  E,  ft'  •  ft,  f)T(r,  ft',  E' ,  t) 


+  qe(f ,  ft,  E,  t)  -  - 


1/d  T{t)V{r,n,E,t)  d 


v\  it  T(t)  (2.83) 


The  advantage  of  decomposing  the  angular  flux  distribution  in  this  manner  is  that 
TDTORT  can  utilize  large  time  steps  when  solving  the  shape  function  and  shorter 
time  steps  when  solving  the  amplitude  function.  I  did  not  implement  this  method  in 
my  code. 


2.5  Time- Dependent  Distribution  Iteration 

The  key  to  DI  is  the  separation  of  the  within  cell  transport  from  the  global 
partial  currents  problem.  One  step  in  affecting  that  separation  is  the  integration  over 
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the  spatial  cell,  which  transforms  the  problem  into  a  coupled  set  of  partial  currents 
that  can  be  solved  exactly.  Starting  from  the  within-group  BTE  in  slab  geometry 

1  d^n(x,t)  &l/}n(x,t)  .  . 

- W7 - h  Hn - X - h  crt(x,  t)ipn(x,  t)  =  qn(x,  t),  (2.84) 

v  at  ox 

where  n  is  the  ordinate  index,  and  averaging  over  the  spatial  cell  yields 


1  fXRld'tpn(x,t)  1  fXR  di/jn(x,t) 

Kxjtl  ^  +  A xJXL 

1  fXR 

+  Ax  J  at(x,t)'iljn(x,t)dx 


1  [XR 

—  J  qn(x,t)dx  ,  (2.85) 


where  Ax  =  xr  —  x l-  Using  the  first  mean  value  theorem  for  integration,  the  cross 
section  term  becomes 


hS,„  <x^x^Ax  =^1* 


i>n{x,t)  dx 


(2.86) 


where  £  G  [xl,xr]  provided  that  a(x,t)  is  continuous  over  the  interval,  which  will  be 
true  if  cells  do  not  span  material  boundaries.  Define  the  cell  average  angular  flux  and 
the  cell  average  source,  q,  as 


1  fXR 

^  /  ^n(x,t)dx 

1  fXR 

9n,i(D  = -x—  /  qn{x,t)dx  , 


(2.87) 


(2.88) 


where  i  indexes  the  cell.  The  BTE  simplihes  to 


1  &lpn(x,t) 


dx  +  XT  “  V£,i(*))  +  VuWnAt)  =  (2-89) 


Ax  JXL  V  dt 


where  ?/;U  and  are  the  angular  flux  distributions  on  the  left  and  right  cell  bound¬ 
aries,  respectively. 
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In  order  to  achieve  a  more  useful  equation,  the  leading  term  needs  the  order  of 
the  integration  and  differentiation  interchanged.  Assuming  that  tfj  meets  the  criteria 
for  the  Lebesgue  dominated  convergence  theorem,  which  requires  that  the  gradient 
is  finite  almost  everywhere  in  the  domain2,  then  (because  v  is  constant  under  the 
multigroup  energy  discretization) 


J_  pl^nMdl 

A  X  JxL  V  & t 


1  d 

vAx  dt 


cxr 


1pn(x,  t)dx  . 


Thus,  substituting  (2.87)  yields 


1  d 

v  dt 


+  ^  WmW  ~  V’n.iW)  +  VtMlpnM 


Qn,i  (t)  • 


(2.90) 


(2.91) 


Applying  the  backward  Euler  method  to  the  above  ordinary  differential  equation, 
yields 


/^7l 

Ax 


(V’n.ij+l  -  *Ph,i,j+ 1)  + 


1  A  ./.A  *  ,1 

ati’j+1  +  vAt )  ^n"+1  -  qn’i’j+1  +  vAt 


n,i,j  ’ 


(2.92) 


where  j  is  the  time  index. 


A  critical  aspect  about  the  synthetic  total  cross  section  approach  is  that  it  makes 
the  transport  problem  more  absorptive  than  it  would  be  otherwise.  The  synthetic 
scattering  ratio  is 


c  = 


&S 


CTb+i 


l 

vAt 


(2.93) 


As  vAt  — >  0,  the  synthetic  scattering  ratio  goes  to  zero.  This  raises  the  question  as 
to  whether  DI  will  maintain  a  performance  advantage  vis-a-vis  source  iteration.  The 
synthetic  total  cross  section  also  makes  the  cells  optically  thick,  which  may  impact 
the  ability  to  determine  the  inward  angular  partial  current  distribution  because  of 
the  attenuation  of  the  angular  partial  currents.  For  short  time  steps,  very  few  neu- 


2Because  this  may  not  be  strictly  true  for  delta  distributed  sources,  a  first  flux  and  first  scatter 
source  approach  may  be  necessary 
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trons  propagate  across  the  cell,  which  makes  the  angular  partial  current  distribution 
susceptible  to  numerical  noise. 


2.5.1  Derivation  of  Diamond  Difference  with  Time  Dependence.  In  this  sec¬ 
tion,  I  will  show  the  derivation  of  the  diamond  difference  spatial  quadrature  transport 
coefficients  with  the  time  derivative  present.  Start  with  the  balance  equation  in  slab 
geometry,  discretized  with  the  backward  Euler  method 


1 

v 


tf+i  ~ 

At 


+  n 


Vf+i  -  rf+i 


Ax 


+  (Ttj+1tf+ 1  =  Qj+ 1- 


(2.94) 


Using  the  auxiliary  equation  to  eliminate  £R  from  the  equation  yields 


•3+1 


+ fk)  u+. + £  (2U«  -  2d«) = 5f+i + (2.95) 


Let  crt  +1  =  crt  +  l/(uAt)  and  ex  =  aAx/fi,  then 


^3+ 1  = 


ex  +  2 


(2.96) 


Substituting  £+i  =  yields 


,A  EQi+if  +  W+i 


^i+i  = 


+  2 


(2.97) 


which  is  the  same  form  for  £A  as  for  the  time-independent  case.  Using  the  synthetic 
total  cross  section  and  synthetic  source,  any  spatial  quadrature  can  be  used  in  a 
time-dependent  transport  code. 


2.5.2  Derivation  of  Step  Characteristic  with  Time  Dependence.  It  is  useful 
to  show  the  derivation  of  a  characteristic  method  in  a  time  dependent  form.  To  derive 
the  step  characteristic  spatial  quadrature  with  the  time  dependence  term  we  begin 
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with  the  BTE  in  slab  geometry  discretized  with  the  backward  Euler  method: 


1  ( '• ’i  1 1  ~l'A  ,  ^+i  i  / 

v  \  At  )  +  +  <T,J+'^+1  ~  «f+‘- 


(2.98) 


Let  (tt  +1  =  cr[j  +  1]  +  l/(r>At)  and  dividing  by  //  yields 


,  ^"b+i 


1  1  Vfi 


ax  +  yr^+1  “  A+1  + 


(2.99) 


As  explained  in  the  previous  section,  the  cross  section  is  constant  within  a  cell,  thus 
at  is  also  constant  within  a  cell.  Define  r(x,x0)  =  (x  —  x0)a tj+1/n  where  x  and  xa  are 
in  the  same  cell.  Multiplying  by  eT(x,x°\  we  have 


<%+ 1 
dx 


er{x,x0)  +  ah+l^.+ier{x,Xo)  =  Agi+ier(*,*o)  +  AAA_eA*,Zo) 


1  A 


H  vAt 


(2.100) 


Without  loss  of  generality,  we  can  consider  the  n  >  0  case.  Noting  that 


d_ 

dx 


(il)er^x'x°>\  —  ^—eT^x,Xo^  +  —il)eT^x'x°> 


(9a; 


(2.101) 


and  integrating  (2.100)  from  x0  =  Xl  to  i  yields 


^+i(x)er^’XL)  -  ^+1(xL)er^L’^  =  [  dx'  (M+lAAeT{x',xA 


'  XL 


+ iy  <^o2> 


Dividing  by  eT(x,XL)  and  rearranging  yields 


^■+i(x)  =  ^■+i(a;L)e-T(3:’a:L)  +  [  dx'  £i±l 


'*L 


+  J  dx'  (2.103) 


Note  that 


t(x',xl)  -t(x,xl)  =  -(x  -  x')atj+1/n 


(2.104) 


41 


can  be  expressed  as 


r(x',xL)  -  t(x,xl)  =  -ix(x  -  x')/Ax, 


(2.105) 


where  ix  =  at.+1Ax//x.  Let  x  =  xr,  then  the  transport  equation  for  a  generic  cell  is 


fpf+i  =  V'j+ie  +  /  dx 


r%  r 


/  Qj+ljx')  „-ex(XTI-x')/Ax 


’XL 


j-X  R 

+  /  dx 

j  #L 


'  1  ^V)e“e‘e(a:R_i'')/Aie 


IxvAt 


.  (2.106) 


For  step  characteristic,  the  underlying  assumption  is  that  the  scattering  and  emission 
source  is  constant  in  a  cell,  thus 


qj+i(x)  =  qj+ 1. 


(2.107) 


Define  the  exponential  moment  function  as  (Mathews  [25]) 


Mn(e)  =  /  (1  -  u)ne~tudu  . 

Jo 


(2.108) 


Using  the  transformation  u  =  (%  —  x') /Ax,  the  source  term  can  be  rewritten  as 


(1^J  [  R  dx'  e-^n~x')/Ax  =  q  ^  [  d u  e 

h  Jx  L  I1  Jo 

_  Ax 

=  qj+ 1 - M  o{ex). 

/i 


(2.109) 

(2.110) 


Note  that  qj+iAx//i  is  the  integrated  production  rate  density  of  particles  through  the 
cell  along  the  ray  and  _Ado(ea;)  is  the  fraction  that  reaches  the  right  boundary.  Making 
the  above  substitution  yields 


Ai .  r* .  ,  i 


V'j+i  =  V’i+ie  €*  +  qj+i—M0(ex)  +  dx ' 

h  (IV  LAI 


(2.111) 
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Using  the  step  approximation  ^ j{x )  =  to  perform  the  following  integration 


'  Ax'  V’j(a:,)e_gx(SR-!e,)/Aa 


(2.112) 


yields 


vAt  /I  JXl  Ax 


vAt  J  ji 


(2.113) 


Thus,  we  have 


A,j+i  =  A,j+ie  €x  +  Um+i  +  ^tA,j  j  —Mo(ex). 


(2.114) 


This  is  the  same  solution  obtained  when  using  the  synthetic  cross  section  and  synthetic 


source. 


The  cell  average  flux  is  defined  as 


rf+i  =  / 


(2.115) 


Substituting  in  (2.103)  and  using  the  fi  >  0  case  produces 


A  _  /  e-ex(x-XL,)/Ax  ,  f  dx^<lj+ 1  f  ,  /  -ex(x-x')/Ax 

* i+1  ~  L  ^i+'  +  XL  Ax  „  XL 


+  [ " A-  f  ix'  (2,ii6) 


Ax  ,7X,  /it’Af 


Using  the  exponential  moment  functions,  this  simplifies  to 


V’j+i  =  ^j+i-Mofe)  H - -gj+i.Mi(ex) 


+  r^.  rdX'Mfie-^-,^,  (2,ii7) 


Ax  ,7X,  /it’Af 
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Next,  the  following  integration  needs  to  be  performed 


I 


c-ex(x-x')/Ax 

/ivA  t 


Using  the  step  approximation  to  solve  the  integral  /  yields 


/\  qp 

^ t+i  =  rf+1M0(ex)  +  —  (  qj+ 1  + 

I-1- 


Aii(ex). 


(2.118) 


(2.119) 


This  is,  again,  the  same  equation  obtained  by  using  the  synthetic  cross  section  and 
synthetic  source. 
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III.  Matrix  Albedo 


While  developing  the  verification  plan  (Chapter  V)  there  was  a  need  for  a 
non-discrete  ordinates  benchmark  that  was  more  accurate  than  a  Monte  Carlo  code. 
Ganapol  [13]  developed  the  TIEL  benchmark,  which  uses  a  Fourier  transform  solution 
based  on  an  analytical  moments  representation  of  the  Green’s  function.  The  bench¬ 
mark  features  a  delta-distributed  source  in  an  infinite  medium.  The  TIEL  benchmark 
is  a  semi-analytic  solution  to 

=  +  (3.1) 

z  1=0 

where  Ui  are  the  scattering  moments,  Pi  are  the  Legendre  polynomials  and  are 
the  Legendre  moments  given  by 

1  f1 

Mx)  =  2  J  d/i  Pj(/j)#e,/j).  (3.2) 

This  benchmark  presents  two  challenging  problems.  First,  the  delta-distributed  source 
must  be  adequately  represented  in  the  transport  code  and,  second,  the  infinite  medium 
has  to  be  properly  handled.  The  other  motivation  for  implementing  this  benchmark 
is  that  it  is  similar  to  the  monoenergetic,  time  dependent  benchmark  also  developed 
by  Ganapol  [12]. 

The  delta-distributed  source  can  be  implemented  in  two  different  fashions.  The 
first  method  would  be  to  introduce  a  single  cell  source  region;  however,  this  is  a  crude 
approach  because  it  effectively  distributes  the  source  into  a  small  volume.  An  alter¬ 
native  approach  would  be  to  use  a  boundary  current  in  conjunction  with  a  symmetry 
boundary  condition. 

Usually,  an  incident  current  is  specified  as  a  Lambertian  source,  i.e.,  the  incident 
flux  is  independent  of  direction.  Alternatively,  one  can  specify  a  source  such  that 
the  incident  current  is  isotropically  distributed  (Appendix  C).  The  delta-distributed 
source  in  the  TIEL  benchmark  is  equivalent  to  an  isotropically  distributed  incident 


a  , 

fix - 1-  1 

ox 
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current.  Isotropic  currents  are  particularly  challenging  because  the  incident  flux  is 
singular  as  /i  — >  0.  The  accuracy  of  the  incident  current  representation,  therefore,  is 
dependent  upon  the  order  of  the  angular  quadrature. 

The  implementation  of  an  incident  current  with  a  non-vacuum  boundary  con¬ 
dition  is  not  a  common  feature  found  in  production  transport  codes.  I  did  implement 
it  in  my  DI  code  because  it  provided  two  different  methods  for  utilizing  the  TIEL 
benchmark. 

3. 1  Infinite  Medium 

While  the  treatment  of  the  delta-distributed  source  is  relatively  straightforward, 
the  infinite  medium  is  not  as  easy.  There  are  traditionally  two  different  approaches 
to  handling  an  infinite  medium.  The  first  method  is  to  define  a  region  much  larger 
than  the  region  of  interest  and  use  vacuum  boundaries.  For  example,  if  the  region 
of  interest  is  4  mean  free  paths  from  the  origin,  the  problem  might  be  defined  as  a 
100  mean  free  path  region.  The  second  method  is  to  utilize  either  a  specular  or  grey 
boundary  condition  and  define  a  region  smaller  than  the  first  method.  Both  of  these 
methods  are  deficient  when  the  scattering  ratio  is  non-zero.  The  first  method  will 
underestimate  the  scalar  flux  because  the  contribution  from  the  medium  beyond  the 
defined  region  is  omitted.  The  reflected  angular  flux  distribution  produced  by  either 
the  specular  or  the  grey  boundary  conditions  will  be  inaccurate  (Figure  3.1). 

None  of  these  methods  were  acceptable  in  the  verification  testing  because  of  the 
difficulty  of  attributing  the  cause  of  any  discrepancies  between  the  results  produced  by 
DI  and  the  TIEL  solution.  Instead,  a  boundary  condition  that  produces  an  accurate 
reflected  flux  is  required  to  reliably  conduct  the  verification  testing.  Such  a  boundary 
condition  could  be  implemented  as  a  matrix,  which  properly  distributes  the  outgoing 
flux  to  match  the  returning  flux  from  an  infinite  medium. 
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Figure  3.1:  Conceptual  flux  distributions  at  a  bound¬ 
ary. 

3.2  Matrix  Albedo  in  Slab  Geometry 

Consider  a  region  within  the  infinite  medium  that  encompasses  any  points  of 
interest.  Without  loss  of  generality,  consider  only  the  right  boundary  by  adding  one 
cell  to  the  right  boundary  of  the  defined  region.  Define  the  angular  flux  vectors 
entering  and  exiting  the  additional  cell  as  shown  in  Figure  3.2.  The  outgoing  fluxes 


Figure  3.2:  Region  of  interest  with  the  one  cell  exten¬ 
sion. 

("0b  and  0^)  are  related  to  the  incoming  fluxes  (0^  and  -0^)  through  the  expression 
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where  Mlr  is  the  transmission  coefficient  matrix  from  the  right  face  to  the  left  face, 
Mll  is  the  reflection  coefficient  matrix  for  the  left  face,  Mrr  is  the  reflection  coefficient 
matrix  for  the  right  face  and  Mrl  is  the  transmission  coefficient  matrix  from  the 
left  face  to  the  right  face.  Recall  that  (2.35)  expresses  the  relationship  between  the 
outgoing  currents  (or  fluxes)  and  incoming  currents  (or  fluxes).  Thus,  the  coefficient 
matrices  are  extracted  from  the  m01  matrix  on  a  flow  direction  basis.  The  columns  of 
the  M  matrices  correspond  to  the  inward  flow  and  the  rows  correspond  to  the  outward 
flow.  The  Mrl  matrix,  for  example,  is  constructed  by  populating  the  elements  where 
the  columns  map  to  the  positive  ordinates  and  the  rows  map  to  the  negative  ordinates. 

Define  a  matrix  albedo  at  the  boundar  of  the  region  of  interest,  A,  such  that 

Vt  =  Hi-  (3-4) 

Because  the  extra  cell  is  one  cell  of  an  infinite  medium,  the  medium  to  the  right  of  the 
right  boundary  of  the  cell  is  identical  to  the  medium  to  the  right  of  the  left  boundary 
of  the  cell,  thus 

VIr  =  A^.  (3-5) 

If  the  incident  flux  is  replaced  with  the  identity  matrix,  which  effectively  is  a 
matrix  constructed  from  the  unit  incident  flux  vectors  for  each  ordinate,  then 


\hL  =  AI  =  A, 


(3.6) 


where  is  the  return  flux  matrix.  Let  F  be  the  incident  flux  matrix  on  the  right 
boundary,  then 

=  AF.  (3.7) 


Substituting  into  (3.3)  yields 


A 

Mlr 

Mll 

AF 

F 

Mrr 

Mrl 

I 

(3.8) 
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One  obvious  method  for  solving  the  above  equation  for  A  is  a  form  of  source  iteration. 
Consider  a  flight  of  neutrons  that  enter  the  cell  from  the  left.  Some  are  reflected  back 
(A0  =  MllI)  and  some  propagate  to  the  right  boundary  (F0  =  MrjJ)  .  A  second 
flight  of  neutrons  enter  the  cell  from  the  left.  Again,  some  neutrons  are  reflected 
plus  there  is  a  contribution  from  the  neutrons  that  reflected  after  exiting  the  right 
boundary  (Ai  =  MllI  +  MlrAoFq).  The  algorithm  that  implements  this  iteration  is 
shown  in  Algorithm  3.1.  An  alternative  iteration  strategy  is  shown  in  Algorithm  3.2. 
I  implemented  both  methods  and  there  was  no  difference  in  performance. 


Compute  U  =  MllMrl 

repeat 

Compute  A  =  MlrU  +  Mll 
Compute  F  =  MrrU  +  MRL 
Compute  U  =  AF 
until  U  converges 

Algorithm  3.1:  Iteration  algorithm  for  determining  the 
matrix  albedo 


Compute  U  =  MllMrl 

repeat 

Compute  A  =  MlrU  +  Mll 
Solve  (I  -  AMrr)U  =  AMrl  for  U 
until  U  converges 

Algorithm  3.2:  Alternative  iteration  algorithm  for  de¬ 
termining  the  matrix  albedo 
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IV.  Implementation 


I  pursued  my  research  effort  in  two  parallel  tracks.  The  first  track  was  to 
develop  a  robust  time  independent  distribution  iteration  code.  The  second  track 
was  to  adapt  the  distribution  iteration  method  into  a  time  dependent  scheme.  The 
essential  element  of  the  first  track  was  the  application  of  modern  software  engineering 
methods  to  a  nuclear  engineering  problem.  The  second  track  required  the  assessment 
of  the  various  mathematical  treatments  of  the  discretization  of  the  time  derivative. 
These  established  the  requirements  for  my  implementation  efforts. 

4-1  Development  of  Time  Independent  Distribution  Iteration  Code 

The  implementation  of  the  time-independent  one-dimensional  slab  and  two- 
dimensional  XY  geometries  followed  the  work  of  Wager  [32]  and  Prins  [28]  closely. 
Prins  presented  a  block  vector  and  matrix  notation  that  captured  the  matrices  and 
vectors  associated  with  X  and  Y  dimensions  into  one  set  of  equations.  I  followed 
a  logical  approach  and  extended  the  notation  to  include  the  vectors  and  matrices 
associated  with  the  Z  dimension. 

4-1.1  Software  Engineering  Practices.  I  quickly  decided  to  start  with  a 
blank  sheet  and  not  reuse  the  software  developed  by  Wager  and  Prins.  The  primary 
reason  was  that  during  the  course  of  their  work,  many  avenues  were  explored  which 
contributed  to  needless  complexity  and  residual  constraints  of  their  implementations 
of  the  DI  method.  Cleaning  the  existing  code  would  have  been  harder  than  starting 
over.  Another  reason  for  starting  over  was  the  goal  to  extend  DI  into  three  dimensional 
XYZ  geometry.  Having  XYZ  geometry  as  a  design  goal  impacts  the  design  of  data 
structures.  The  third  reason  was  the  fact  that  a  time  dependent  algorithm  was  going 
to  be  implemented  using  DI. 

In  the  development  of  my  DI  code,  I  adopted  several  key  practices  that  are  used 
in  the  software  engineering  field,  including: 

1.  Using  a  source  code  management  system; 
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2.  Testing  to  ensure  that  revisions  do  not  prevent  compilation  or  introduce  error; 

3.  Committing  changes  to  the  source  management  system  frequently. 

The  two  key  features  that  I  wanted  from  a  source  code  management  system  were  the 
ability  to  branch  and  to  rollback  changes.  The  ability  to  branch,  that  is  to  have  a 
separate  development  effort,  from  the  main  development  effort  is  an  important  feature 
to  have  in  order  to  support  concurrent  research  efforts.  The  second  and  third  practices 
resulted  in  a  coding  work  flow  that  focused  on  making  small  incremental  changes  and 
frequent  testing. 

Based  on  the  development  practices  that  I  adopted,  the  strategy  for  developing 
the  DI  code  was  to  first  implement  a  slab  geometry  version  and  verify  its  performance. 
Once  a  working  implementation  of  slab  geometry  was  available,  even  though  the 
feature  set  was  incomplete,  the  functional  slab  geometry  was  branched  from  the  main 
development  effort.  The  code  was  then  reviewed  for  the  elements  common  between 
XY  and  slab  geometry.  The  common  elements  were  separated  from  the  geometry 
specific  elements  and  two  libraries  were  created.  An  XY  geometry  specific  library 
was  created,  using  the  slab  geometry  specific  library  as  a  template.  This  resulted 
in  a  functional  code  that  consisted  of  a  verified  slab  geometry  component  and  an 
unverified  XY  geometry  component. 

The  main  trunk  development  at  this  point  was  focused  on  verifying  the  perfor¬ 
mance  of  the  XY  geometry  and  correcting  errors.  The  test  suite  consisted  of  both 
slab  and  XY  geometry  cases  and  changes  were  reviewed  to  make  sure  slab  geome¬ 
try  performance  was  not  affected.  Additional  features  were  added  to  slab  geometry 
to  support  research  efforts  outside  the  scope  of  this  dissertation.  Support  for  time 
dependent  transport  was  added  and  testing  in  slab  geometry  was  performed.  The 
introduction  of  time  dependent  support  at  this  juncture  was  done  to  ensure  that  the 
evolution  of  the  source  code  would  be  compatible  with  the  goals  of  my  research. 

Once  a  functional  and  verified  XY  geometry  implementation  was  available,  the 
slab  and  XY  geometry  implementation  was  branched  and  development  of  a  XYZ  ge- 
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ometry  implementation  was  started  in  the  main  trunk.  The  XY  geometry  specific 
library  was  used  as  a  template  for  the  XYZ  geometry.  During  the  course  of  the 
development  of  the  XYZ  geometry  version,  additional  code  elements  that  would  be 
common  to  all  three  geometries  became  evident;  however,  1  did  not  attempt  to  ex¬ 
tract  these  newly  identified  common  elements  because  they  required  some  substantial 
changes  to  the  code  base.  In  keeping  with  the  development  principles  that  I  adopted, 
the  extraction  of  the  common  elements  would  require  a  parallel  development  effort 
and  then  a  subsequent  merge  of  the  changes  into  the  main  development  trunk.  The 
extraction  of  these  common  elements  was  deferred  at  this  time;  however,  I  recom¬ 
mend  that  it  be  pursued  prior  to  any  work  in  developing  a  message  passing  interface 
parallelized  version. 

4-2  Three  Dimensional  Geometry 

The  first  step  in  developing  the  three-dimensional  implementation  of  DI,  or  any 
transport  code  for  that  matter,  is  to  establish  the  cell  orientation.  Because  1  wanted  to 
maintain  the  orientation  where  the  XY  plane  is  normal  to  an  observer’s  line  of  sight, 
the  Z  axis  extends  towards  the  observer  in  the  positive  Z  direction.  This  arrangement 
forms  the  basis  of  the  arrangement  of  the  cell  faces,  as  shown  in  Figure  4.1.  While 
this  arrangement  is  purely  arbitrary,  it  does  facilitate  code  reuse  with  the  XY  and 
slab  geometries. 


Figure  4.1:  Three-Dimensional  Cube  Orientation  and 
Face  Naming 

Given  the  arrangement  of  faces  of  a  computational  cell,  the  next  step  is  to 
establish  an  ordering  convention  for  use  when  constructing  the  partial  currents  linear 


52 


system.  Again,  I  kept  the  order  established  in  the  XY  geometry  and  added  the  Z 
faces  after  the  XY  faces.  This  results  in  an  arrangement  where  all  the  X  faces  (YZ 
planes)  are  first,  then  the  Y  faces  (XZ  planes)  and  finally  the  Z  faces  (XY  planes). 
When  implementing  DI,  one  has  to  decide  on  whether  to  use  angular  fluxes  or  angular 
currents  to  represent  flows  into  and  out  of  the  cell.  I  used  the  current  representation 
because  Prins  [28]  found  that  DI  converges  more  rapidly  with  it  than  with  the  angular 
flux  representation. 


4-2.1  Zeroth  Spatial  Moment  Methods.  In  XYZ  geometry,  I  adopted  the 
block  vector  and  matrix  notation  used  by  Prins  [28].  The  cell  face  angular  current 
vector  is 


J  =  Ul  Jr  3b  3t  Jp  3f ) 


(4.1) 


where  jp  is  the  current  on  the  left  face,  jp  is  the  current  on  the  right  face,  jp  is  the 
current  on  the  bottom  face,  jp  is  the  current  on  the  top  face,  jp  is  the  current  on  the 
back  face1  and  jp  is  the  current  on  the  front  face.  The  scattering  and  emission  source 
in  a  cell  for  a  zeroth  moment  method  is  simply 


qs  =  f 


(4.2) 


and 

qe  =  r  (4.3) 

where  (f  is  the  cell  average  scattering  source  and  qe  is  the  cell  average  emission  source. 
The  system  of  discrete  ordinates  spatial  quadrature  equations  for  within-cell  transport 
is 

Jout  =  KOIJm  +  KOSgs  +  KOEge  (4.4) 

and 

_  Kqijin  +  K^Sgs  +  (4.5) 

1The  superscript  “P”  is  from  “Posterior”  or,  if  you  prefer,  the  Greek  translation  of  back  is  ttujuj. 
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As  in  the  XY  geometry  case,  we  define  the  augmented  K01  matrix: 

K™  K“  K™  K“  K“ 

^BL  ^ RR  ®^BB  ^RT  ^RP  ^RP 

^BL  ®^BR  ^BT  ^BP  ®BF 

K“  K^b  K“  K^p 

Kp\  K™  K%  K%  K% 

K^r  K fT  K% 

The  subscripts  denote  the  combination  of  cell  faces  using  the  notation  (to,  from),  e.g., 
the  subscript  “LT':  is  for  flows  to  the  Left  face  from  the  Top  face.  The  submatrices, 
e.g.  1K(^,  operate  on  a  subvector  of  angular  partial  currents  and  are  sized  by  the 
number  of  ordinates.  The  K01  matrix  has  a  very  sparse  structure  because  there  is 
no  coupling  between  ordinates  as  a  result  of  scattering.  All  the  submatrices  along 
the  main  diagonal,  e.g.  are  zero  matrices.  The  off-diagonal  submatrices,  e.g. 

only  have  non-zero  elements  where  there  is  a  coupling  between  an  ordinate  on 
the  inflow  face  and  an  ordinate  on  the  outflow  face.  The  location  of  the  non-zero 
elements,  as  well  as  the  total  number,  in  each  submatrix  depends  upon  the  angular 
quadrature  used.  For  example,  in  the  Sn  angular  quadrature  there  are  n(n  +  2) 
ordinates,  which  results  in  a  IK01  matrix  with  dimensions  of  6 n(n  +  2)  x  6 n(n  +  2). 
Each  off-diagonal  submatrix  has  n(n  +  2)/2  non- zero  elements  along  its  main  diagonal 
for  a  total  of  15n(n  +  2). 

The  Kos  and  IKOE  matrices  are 

K°s  _  (Kos  Kos  Kos  Kos  Kos  KosjT  (4.7) 

Koe=(k(Ie  KgE  KgE  K?e  K%e  K°eY  ■  (4.8) 

For  the  above  matrices,  the  subscript  denotes  the  outflow  face.  Both  the  Kos  and  KOE 
matrices  are  sparse-the  non-zero  elements  are  where  an  ordinate  is  in  the  outward 
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direction  for  a  face.  For  example,  the  IK submatrix  has  non-zero  elements  for 
ordinates  where  fi  <  0. 

The  IK^1  matrix  is  defined  as 


.  (4.9) 


In  the  case  of  the  IK^1  matrix,  the  subscript  denotes  the  inflow  face.  The  KA1  matrix  is 
sparse  as  each  submatrix  is  diagonal.  For  zeroth  moment  methods,  the  KW  and 
matrices  are  the  same  in  slab,  XY  and  XYZ  geometries.  The  augmented  scattering 
matrix  is  defined  as 


/ 


X  = 


Xs  0  0 

0  Xs  0 
0  0  X, 


\ 


) 


(4.10) 


where  Xs  is  the  isotropic  scattering  matrix.  For  anisotropic  scatter  each  submatrix 
will  be  different.  The  Xs  matrix  is  not  sparse,  though  the  augmented  matrix,  X,  is 
sparse. 


While  the  augmented  IK  and  X  matrices  are  large,  the  primary  benefit  of  using 
them  is  that  the  implementation  becomes  simpler.  The  large  dimensions  of  the  ma¬ 
trices  can  be  mitigated  by  representing  them  as  sparse  matrices  and  using  a  sparse 
Basic  Linear  Algebra  System  library.  For  example,  in  the  S§  angular  quadrature,  the 
IK01  matrix  is  480  x  480  with  2000  non-zero  elements,  which  is  less  than  one  percent 
of  all  the  elements.  The  IKOE  and  IKos  are  480  x  80  with  240  non-zero  elements  and 
the  KF’1  matrix  is  80  x  480  with  480  non-zero  elements.  The  X  matrix  is  the  most 
dense  matrix  with  a  third  of  its  elements  being  non-zero. 


The  solution  of  the  within-cell  transport  equations  follows  directly  from  Math¬ 
ews  [24],  By  using  the  augmented  vector  and  matrix  notation,  we  have  the  same  set 
of  equations  in  slab,  XY,  and  XYZ  geometry. 
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4-2.2  First  Spatial  Moment  Methods.  All  the  vectors  and  matrices  can  be 
further  expanded  by  adding  the  components  for  the  first  spatial  moments  in  a  block 
fashion  as  demonstrated  by  Prins  [24],  The  use  of  the  block  notation  allows  the 
implementation  of  a  unified  code  that  handles  both  zero  and  first  spatial  moment 
methods. 

4-2.3  Boundary  Conditions.  Boundary  conditions  can  either  be  explicit  or 
implicit.  An  explicit  boundary  condition  provides  a  known  incoming  flux  (or  current) 
on  the  surface  of  the  problem  domain.  An  example  of  a  explicit  boundary  condi¬ 
tion  is  a  Lambertian  illumination  (Appendix  C).  An  implicit  boundary  condition 
redistributes  the  outward  flux  back  into  the  problem  domain. 

An  incident  illumination  can  either  be  specified  as  an  angular  flux,  which  is  the 
typical  approach,  or  as  a  current.  My  D1  code  supports  both  notations  to  facilitate 
my  verification  effort.  Because  1  used  the  current  representation,  angular  fluxes  need 
to  be  converted  into  partial  currents  (both  J  and  j).  For  simplicity  and  to  facilitate 
angular  quadrature  refinement  testing,  1  only  implemented  an  isotropic  angular  flux 
(a  Lambertian  illumination) — an  anisotropic  angular  flux  is  a  trivial  modification  of 
the  input  hie  processor.  For  an  incident  current,  the  partial  current  and  a  distribu¬ 
tion  (Lambertian  or  Isotropic  Surface  Source)  is  specified  as  input.  The  code  then 
computes  the  angular  partial  current  using  the  specified  distribution. 

Given  an  incident  angular  flux  ^(fl),  the  inward  partial  current  on  face,  J™  is 
given  by 

Jn  =  f  do  I h  ■  n\ip(Cl)  (4.11) 

J  h-Ci<  o 

where  the  subscript  n  denotes  the  face  and  n  is  the  outward  unit  vector  normal  to 
the  face.  The  inward  current  on  a  cell  is  j^n(0)  =  | n  ■  fl|^(fl). 


56 


4-3  Implementing  the  Matrix  Albedo  Into  Distribution  Iteration 

The  matrix  albedo  was  implemented  in  my  DI  code  for  slab  geometry.  Several 
changes  are  required  to  implement  the  matrix  albedo: 

1.  Calculate  the  transport  coefficients,  e.g.  IK01,  for  the  infinite  medium  material. 

2.  Calculate  the  m01  matrix  for  the  infinite  medium  material. 

3.  Implement  a  solver  for  the  albedo  matrix  A. 

4.  Implement  all  boundary  conditions,  e.g.  specular  reflection,  in  matrix  form  and 
use  them  when  improving  the  inflow  flux  distribution. 

5.  Provide  a  routine  to  compute  the  reflection  coefficient  used  in  the  partial  current 
solver. 

Steps  1  and  2  were  easy  to  implement  in  my  DI  code  because  the  code  necessary 
to  compute  the  coefficient  matrices  and  m01  was  already  implemented.  The  code  was 
implemented  to  permit  the  use  of  different  materials  on  the  left  and  right  boundaries. 
For  example,  a  three  region  problem  could  have  an  infinite  region  of  material  “A”  on 
the  left,  a  region  of  interest  of  finite  length  of  material  “B”  and  an  infinite  region  of 
material  “C”  on  the  left. 

I  implemented  two  different  solvers  for  the  matrix  albedo  (Mathews  &  Dishaw 
[23]).  Because  both  solvers  are  iterative  solvers,  care  must  be  taken  to  avoid  false 
convergence.  Both  algorithms  check  for  convergence  using  an  iteration  count  doubling 
strategy  shown  in  Algorithm  4.1.  This  algorithm  produced  reliable  results  even  in  a 
strongly  scattering  medium,  e.g.  c  =  0.9. 

In  order  to  facilitate  the  use  of  specular  and  grey  boundaries  along  with  the  infi¬ 
nite  medium  matrix  albedo  boundary,  all  boundaries  were  represented  with  a  matrix. 
A  specular  boundary  can  be  represented  in  matrix  form  as  (where  ordinate  1  and  n 
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Set  iteration  count  to  1 
Set  check  convergence  to  2 
Initialize  Up  to  zero 

repeat 

Solve  for  U,  e.g.  Algorithm  3.1 

if  iteration  count  is  equal  to  check  convergence  then 

if  Symmetric  relative  difference  of  U  and  Up  less  than  tolerance  then 
Exit  the  iteration 

end  if 

Up  =  U 

Double  the  value  of  check  convergence 

end  if 

Increment  iteration  count 
until  Maximum  number  of  iterations  are  exceeded 
if  Maximum  number  of  iterations  were  exceeded  then 
Signal  failure  to  converge 

end  if 

Algorithm  4.1:  Iteration  count  algorithm  used  to  de¬ 
termine  if  the  matrix  albedo  has  con¬ 
verged. 


form  a  reflection  pair,  2  and  n  —  1  form  another  reflection  pair,  etc.) 
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where  a  is  the  reflection  coefficient  (a  =  1  is  a  symmetry  boundary,  a  <  1  is  a  dirty 
mirror).  A  grey  boundary  would  have  each  row  of  the  matrix  filled  in. 

The  matrix  albedo  requires  the  implementation  of  a  method  to  determine  the 
reflection  coefficient  that  will  be  used  in  the  partial  current  solver.  The  reflection 
coefficient  at  a  face  is  defined  as 


J 


face 


T  OUt  5 

^  face 


&  face 


(4.13) 


where  Jout  is  the  partial  current  leaving  the  region  and  Jm  is  the  partial  current  that 
has  been  reflected  back.  Let  A face  be  the  matrix  albedo  then,  by  definition, 

Jm  =  w  ■  AfaceZfaceJout,  (4.14) 

where  w  is  a  vector  of  ordinate  weights  and  Z/ace  is  the  angular  current  distribution 
for  a  face,  which  is  a  vector.  The  reflection  coefficient  is,  therefore, 

ol  w  A±faceTjface.  (4.15) 

4-4  Multigroup  and  Time  Dependent  Transport 

The  conventional  approach  for  implementing  time  dependent  and  multigroup 
transport  is  through  an  outer-iteration  method.  For  example,  in  a  time-independent 
multigroup  code  the  outer-iteration  starts  at  the  highest  energy  group  and  uses  a 
time- independent  monoenergetic  transport  solver  to  determine  the  angular  (or  scalar) 
fluxes  for  each  energy  group  using  the  other  group  angular  fluxes  as  part  of  the 
external  source  for  the  within-group  problem.  For  a  problem  where  neutrons  only 
scatter  downwards  in  energy,  only  one  outer-iteration  is  required  because  the  lower 
energy  groups  do  not  couple  to  the  higher  energy  groups.  A  problem  that  includes 
the  upscatter  of  neutrons  or  fission  neutrons,  multiple  outer-iterations  are  required 
because  of  the  coupling  of  the  energy  groups. 

Using  the  synthetic  total  cross-section  and  synthetic  source  method,  time-dependent 
problems  can  be  solved  using  an  outer-iteration  method.  LInlikc  multigroup  trans¬ 
port,  only  one  iteration  for  all  the  time  steps  is  required  because  there  is  no  scat¬ 
tering  of  neutrons  to  previous  time  steps.  A  combined  time-dependent,  multigroup 
transport  code  would  then  consist  of  two  outer-iterations,  with  the  time-dependent 
outer-iteration  calling  the  multigroup  iteration. 

The  implementation  of  a  multigroup  transport  code  involves  a  mechanism  for 
storing  the  group  angular  flux  at  each  time  step,  an  iteration  that  solves  each  energy 
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group  until  all  the  group  angular  fluxes  are  converged,  and  a  subroutine  for  calculating 
the  within-gronp  external  source.  When  using  the  DI  method,  storing  the  inward 
angular  partial  current  distribution  for  each  energy  group  as  the  initialization  (vice 
isotropic)  will  offer  a  performance  advantage. 

The  implicit  backward  Euler  method  serves  as  a  useful  starting  point  for  adapt¬ 
ing  the  Dl  method  into  a  time-dependent  algorithm.  As  shown  in  Section  2.5,  the  Dl 
method  can  be  derived  with  the  time-dependent  form  of  the  Boltzmann  Transport 
Equation.  This  method  does  not  require  any  additional  storage  if  careful  use  of  the 
cell  average  angular  flux  data  structure  is  observed.  The  midpoint  average  method 
that  PARTISN  utilizes  does  require  a  second  copy  of  the  cell  average  angular  flux 
data  structure,  thus  the  price  in  memory  has  to  be  weighed  against  the  improved 
accuracy  in  time.  As  in  the  multigroup  implementation,  a  mechanism  for  storing  and 
retrieving  the  inward  angular  partial  current  distribution  for  the  previous  time  step 
will  offer  a  performance  advantage. 
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V.  Verification 


To  support  my  research,  a  formal  approach  to  verification  of  the  D1  method 
was  needed.  The  primary  reason  for  pursuing  formal  verification  was  to  provide 
a  measure  of  confidence  in  the  underlying  code  base  before  implementing  the  time 
dependent  version  of  DI.  The  secondary  reason  was  to  build  upon  the  work  performed 
by  Wager  [32]  and  Prins  [28]  in  establishing  confidence  in  the  Dl  method  as  a  viable 
alternative  to  the  more  established  source-iteration  based  methods. 

In  the  vernacular  of  the  software  testing  field,  the  formalized  testing  of  software 
is  described  as  “validation  and  verification”  (V&V).  Validation  is  a  measure  of  the 
software’s  ability  to  perform  in  its  intended  application,  e.g.  the  model  is  an  adequate 
analog  to  reality  for  the  intended  application.  Verification  is  a  measurement  of  the 
accuracy  and  reliability  of  the  implementation  of  a  model.  Verification  does  not 
ascertain  whether  or  not  the  implementation  is  suitable  for  an  application.  The 
fidelity  of  the  model  combined  with  the  verification  of  its  implementation  provides  an 
expectation  of  whether  the  combination  is  valid  for  a  specified  application. 

For  the  purpose  of  my  research,  I  focused  on  verification  and  did  not  pursue 
validation.  Because  the  discrete  ordinates  method  which  underlies  DI  has  been  studied 
and  analyzed  for  forty  years  and  its  applicability  is  well  understood  in  the  nuclear 
engineering  community.  Furthermore,  the  focus  of  my  research  was  to  develop  a  more 
robust  implementation  of  the  discrete  ordinates  method  rather  than  to  design  towards 
a  specific  application  of  discrete  ordinates. 

There  is  a  hierarchy  of  models  that  is  intimately  part  of  a  V&V  effort.  First, 
there  is  the  physical  model  which  serves  as  the  description  of  some  physical  process. 
Next  is  the  mathematical  model  that  expresses  the  physical  model.  The  mathematical 
model  can  be  expressed  in  a  form  to  facilitate  its  numerical  solution,  e.g.,  discretiza¬ 
tion.  Finally  there  is  the  computer  model,  which  is  the  model  that  is  ultimately 
implemented  in  software.  Any  one  of  these  models  can  be  the  subject  of  validation 
or  verification. 
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The  purpose  of  verification  is  to  identify  deviations  between  the  output  of  the  im¬ 
plementation  of  a  model  and  the  expected  values  [31].  The  source  of  these  deviations 
can  be  attributed  to  uncertainty  and  error.  Uncertainty  can  be  further  categorized  as 
reducible  (epistemic)  and  irreducible  (aleatoric).  Epistcmic  uncertainty  is  the  result 
of  insufficient  information  about  the  physical  system.  Aleatoric  uncertainty  is  due  to 
the  probabilistic  distribution  of  the  inputs,  e.g.  variation  in  the  material  properties. 
Because  the  goal  of  my  research  is  to  develop  an  improved  algorithm  for  solving  dis¬ 
crete  ordinates,  the  effects  of  epistemic  and  aleatoric  uncertainty  are  irrelevant  and 
need  to  be  eliminated.  Specifically,  my  research  is  focused  on  verifying  that  DI  does 
not  introduce  error  in  the  solution  of  the  discrete  ordinates  equations  for  a  given  spa¬ 
tial  and  angular  quadrature.  In  the  aforementioned  hierarchy  of  models,  my  research 
is  comparing  the  computer  model,  DI,  to  the  numerical  model.  Thus,  all  sources 
of  uncertainty  must  be  eliminated;  for  example,  aleatoric  uncertainty  is  eliminated 
through  the  use  of  specified  material  properties  rather  than  the  use  of  real  material 
properties.  This  permits  me  to  identify  any  sources  of  error  in  my  implementation 
rather  than  attempting  to  discriminate  between  deviations  due  to  uncertainty  and 
those  due  to  error.  Because  DI  is  based  on  discrete  ordinates — thereby  inheriting  its 
benefits  and  deficiencies — the  use  of  a  discrete  ordinates  based  benchmark  is  not  a  sig¬ 
nificant  shortcoming.  If  experimental  data  were  used  instead,  the  effects  of  epistemic 
and  aleatoric  uncertainty  would  need  to  be  identified. 

The  verification  approach  that  I  implemented  is  a  comparison  with  another 
implementation  of  the  discrete  ordinates  equations,  specifically  PARTISN  version  4.00 
(beta  release  05-26-04) 1  (Alcouffe  [2]).  A  suite  of  test  problems  was  constructed  and 
a  suitable  benchmark  was  selected.  While  the  preferred  benchmark  is  an  analytic 
solution,  there  is,  unfortunately,  a  paucity  of  analytic  benchmarks  for  the  Boltzmann 
Transport  Equation.  I  categorize  the  benchmarks  into  the  following  classes: 

•  Analytic  solutions, 

1As  of  21  July  2007  this  was  the  only  version  available  for  distribution.  Los  Alamos  National 
Laboratory  has  not  made  a  newer  version  available 
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•  Semi-analytic  solution  (e.g.  Ganapol’s  benchmarks  [12]  and  [13]), 

•  Exact  discretized  solutions, 

•  Problems  constructed  to  have  specified  solutions,  i.e.,  manufactured  problems, 
and 

•  Established  codes  (e.g.  PARTISN). 

The  use  of  analytic  and  semi-analytic  benchmarks  is  not  necessarily  required  for 
the  type  of  verification  that  I  performed  for  my  research.  The  analytic  and  semi- 
analytic  benchmarks  are  for  verification  of  the  underlying  numerical  model-discrete 
ordinates-to  the  mathematical  model-the  Boltzmann  Transport  Equation.  I  included 
them  because  I  wanted  to  confirm  that  1  was  using  PARTISN  correctly  and  was 
not  introducing  a  common  error.  I  did  not  use  any  manufactured  problems  in  the 
verification  of  my  DI  code. 

5.1  Time- Independent  Verification  Tests 

The  suite  of  tests  is  first  divided  into  the  three  supported  geometries  (Slab, 
2D  rectangular,  and  3D  boxoid).  Grouping  by  geometry  is  required  as  the  different 
types  of  geometries  introduce  different  testing  requirements.  The  tests  are  further 
subdivided  as  to  whether  the  intent  is  to  test  a  qualitative  or  quantitative  aspect.  This 
subdivision  is  more  for  software  development  convenience  than  any  other  reason.  The 
qualitative  tests  were  designed  with  automation  in  mind  so  that  they  could  be  run 
frequently,  typically  before  the  source  code  was  checked  into  the  revision  management 
system  and  also  during  automated  testing  of  a  neutral  build.  This  allowed  for  a 
quick  determination  of  whether  the  software  development  effort  was  on  track.  The 
quantitative  tests  were  designed  to  assess  the  correctness  of  the  DI  code.  These  tests 
required  greater  analytical  effort  and  were  not  fully  automated. 

Wager  [32]  and  Prins  [28]  both  demonstrated  the  performance  of  the  Distribu¬ 
tion  Iteration  method  for  a  series  of  time-independent  problems.  Wager  developed 
DI  originally  and  examined  its  performance  in  slab  geometry  and  Prins  extended  DI 
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to  XY  geometry.  Both  authors  utilized  a  combination  of  analytic  solutions  and  a 
traditional  source  iteration  code  as  benchmarks.  The  suite  of  tests  was  extended  be¬ 
yond  what  was  required  for  verification  in  order  to  evaluate  performance  and  develop 
confidence  in  DI. 

5.1.1  Slab  Geometry.  The  performance  of  the  time  independent,  slab  geom¬ 
etry  DI  method  was  evaluated  by  Wager  [32],  Wager  examined  two  test  cases,  the  first 
one  was  a  single  region  problem  and  the  second  one  was  a  periodic  interface  problem 
consisting  of  two  regions  repeated  ten  times.  Building  from  his  work,  I  implemented 
a  series  of  tests  designed  to  detect  common  programming  errors  and  implementation 
failures.  The  set  of  test  problems  is  presented  in  the  following  list-the  details  of  which 
follow  after  the  list  (dimensions  are  in  units  of  mean  free  path  (MFP)): 

1.  Pure  Absorber-Single  region  (<x  =  l,c  =  0,ge  =  4.7,  32  mfp)  and  vacuum 
boundaries; 

2.  Pure  Scatterer-Single  region  (a  =  l,c  =  l,qe  =  4.7,  32  mfp)  and  vacuum 
boundaries; 

3.  Infinite  Medium-Single  region  (a  =  l,c  =  0.5,  qe  =  4.7,  32  mfp),  symmetry 
boundaries; 

4.  Constant  Source-Single  region  (a  =  0.1,  c  =  0.5,  qe  =  4.7,  3.2  mfp)  vacuum 
boundaries; 

5.  Left  Isotropic  Surface  Source  ( Q  =  0.5),  4  mfp-Single  region  (cr  =  1,  c  =  0.9,  qe  = 
0,  4  mfp),  left  symmetry  boundary  and  right  vacuum  boundary; 

6.  Left  Isotropic  Surface  Source  ( Q  =  0.5),  Infinite  Medium-Single  region  (a  = 
1,  c  —  0.9,  qe  =  0,  Infinite),  left  symmetry  boundary  and  right  albedo  matrix; 

7.  Left  Lambertian  Illumination  ( Q  =  4.7),  4  mfp-Single  Region  (a  —  l,c  — 
0.9,  qe  =  0,  4  mfp),  vacuum  boundaries; 
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8.  ID  Periodic  Interface,  Cross-sections:  Region  A  (cr  =  10,  c  =  1.0,  <7°  =  0,  10 
mfp)  Region  B  (a  =  0.1,  c  =  1.0,  qe  =  0,  0.1  mfp),  Left  Lambertian  illumination 
(Q  =  4.7),  vacuum  boundaries; 

9.  ID  Periodic  Interface,  Scattering  ratio:  Region  A  (a  =  1,  c  =  1,  qe  =  0,  1  mfp) 
Region  B  (a  =  1,  c  =  0,  qe  =  0,  1  mfp),  Left  Lambertian  illumination  (Q  =  4.7), 
vacuum  boundaries; 

10.  ID,  Two  Region,  Scattering  ratio  test,  Region  A  and  B  (a  —  1.0,  S  —  0,  2  mfp), 
vacuum  boundaries,  left  incident  Lambertian  illumination. 

These  benchmark  problems  form  the  basis  of  the  different  test  scenarios  as  shown  in 
Table  5.1.  In  all  cases,  unless  otherwise  noted,  the  convergence  criterion  was  that  the 
symmetric  relative  difference  in  (  between  iterations  was  less  than  lx  ICR7  everywhere. 
The  symmetric  relative  difference  is  defined  as 

SRD(i,s/)  =  2;h^fL.  (5.1) 

M  +  \y\ 

The  key  verification  metric  that  I  evaluated  was  the  maximum  symmetric  relative 
difference  in  the  scalar  flux  computed  by  DI  and  the  benchmark.  If  the  maximum 
symmetric  relative  difference  between  the  two  solutions  was  less  than  ICR6,  then  the 
two  solutions  were  in  agreement.  I  also  examined  the  rate  of  convergence  as  the  spatial 
mesh  was  refined  for  one  of  the  test  problems. 

Because  much  of  this  work  was  performed  by  Wager,  the  focus  of  my  effort  was 
to  verify  that  I  had  correctly  implemented  the  DI  algorithm  rather  than  a  detailed 
assessment  of  the  algorithm.  I  did  expand  upon  Wager’s  assessment  of  DPs  perfor¬ 
mance  by  comparing  the  performance  of  DI  with  a  convergence-accelerated  source 
iteration  code. 

5. 1.1.1  Tests  1  and  2.  The  single  region  pure  absorber  test  (1)  is  the 
simplest  possible  test  and  provides  a  useful  starting  point  for  developing  a  transport 
code.  The  key  advantage  of  this  problem  is  the  fact  that  it  has  an  analytic  solution, 
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Table  5.1:  Summary  of  time  independent,  slab  geometry 
tests. 


Test  Problem 

Analytic 

Exact 

PARTISN 

Pure  Absorber 

All 

DD/LD 

Pure  Scatterer 

DD/LD 

Infinite  Medium 

All 

DD/LD 

Constant  Source,  Single-Region 

SC/LC/LD 

DD/LD 

Isotropic  Source,  4  mfp 

SC 

Isotropic  Source,  Infinite  Medium 
Lambertian  Source,  4  mfp,  Vacuum 

sc 

DD/LD 

Periodic  Interface,  Cross-sections 

DD/LD 

Periodic  Interface,  Scattering  ratio 

DD/LD 

which  greatly  facilitates  the  debugging  of  a  transport  code.  The  complement  to  the 
pure  absorber  test  is  the  pure  scatterer  test  (2).  The  primary  utility  of  the  pure 
scatterer  test  is  that  it  tests  a  different  part  of  the  transport  code,  which  is  useful  for 
debugging  purposes.  Using  a  spatial  mesh  of  256  uniformly  spaced  cells  and  the  Ss 
and  DP4  angular  quadratures  the  solution  produced  by  DI  in  both  tests  agreed  with 
the  benchmark  solutions. 

I  wanted  to  determine  if  there  was  a  difference  in  accuracy  between  the  single- 
range  and  double-range  angular  quadratures.  Carlson  [8]  asserts  that  the  double-range 
angular  quadrature  yields  more  accurate  results  for  thin  cells  than  the  single-range 
angular  quadrature.  Conversely,  the  single-range  angular  quadrature  yields  more  ac¬ 
curate  results  for  thick  cells  because  there  is  a  higher  probability  of  neutrons  traveling 
out  of  a  cell  in  directions  near  /x  =  ±1  and  the  single-range  angular  quadrature  has 
ordinates  closer  to  /i  =  ±1  than  has  the  double-range  angular  quadrature. 

To  determine  if  there  was  such  a  difference  between  the  two  angular  quadratures, 
I  used  a  variant  of  the  pure  absorber  test.  Instead  of  a  uniform  source,  I  used  a 
Lambertian  illumination  on  the  left  boundary.  The  scalar  flux  on  a  cell  face  for  a 
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Lambertian  illumination  is  simply 


<t>(x)  =  j  (e  I<T  -  xar(0,  xa))  ,  (5.2) 

where  r(0,:r<7)  is  the  upper  incomplete  Gamma  function.  LTsing  the  Sg  and  DP4  an¬ 
gular  quadratures,  I  evaluated  the  symmetric  relative  difference  between  the  resulting 
cell  face  scalar  fluxes  and  the  exact  solution.  When  using  step  characteristic,  the 
single-range  angular  quadrature  did  perform  better  when  the  cell  size  was  greater 
than  3  mean  free  paths  (Figure  5.1).  The  linear  discontinuous  spatial  quadrature, 
which  is  not  a  positive  spatial  quadrature,  did  not  exhibit  as  large  of  a  performance 
difference  between  the  two  angular  quadratures.  In  fact,  both  angular  quadratures 
were  ineffective  when  the  cell  size  was  large  enough  to  cause  negative  fluxes  (>  2 
mfp).  For  both  spatial  quadratures,  the  double-range  angular  quadrature  did  have 
superior  accuracy  for  the  thinner  spatial  cells,  as  asserted  by  Carlson. 

5. 1.1. 2  Test  3.  The  infinite  medium  test  (3)  serves  as  an  intermediate 
test  between  the  pure  absorber  and  pure  scatterer  tests.  The  use  of  left  and  right 
symmetry  boundaries,  as  well  as  having  c  ^  1,  yields  an  analytic  solution  for  the  scalar 
flux:  0  =  S/aa.  Using  a  spatial  mesh  of  256  uniformly  spaced  cells  and  the  Sg  and 
DP4  angular  quadratures  the  solutions  produced  by  DI  agreed  with  the  benchmark 
solutions. 


5. 1.1. 3  Test  4-  The  constant  source,  single  region  problem  (4)  served 
as  the  basis  for  evaluating  the  convergence  rate.  Moderate  scattering  (c  =  0.5)  was 
chosen  to  permit  rapid  convergence  of  unaccelerated  source  iteration.  The  spatial 
mesh  was  refined  from  4  cells  to  65536  cells  and  the  test  was  performed  using  the 
various  spatial  quadratures  (SC,  DD,  LD,  and  LC).  Both  single-range  and  double- 
range  Gauss-Legendre  angular  quadratures  were  used  with  8  ordinates  (Sg  and  DP4) 
and  the  convergence  tolerance  was  set  to  10-7.  Using  a  spatial  mesh  of  256  uniformly 
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Figure  5.1:  Comparison  of  the  angular  quadratures  for 
the  step  characteristic  and  linear  discontin¬ 
uous  spatial  quadratures. 
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spaced  cells  and  the  S8  and  DP4  angular  quadratures  the  solutions  produced  by  DI 
agreed  with  the  benchmark  solutions. 

The  convergence  rate  for  each  spatial  and  angular  quadrature  combination  was 
determined  by  computing  the  maximum  symmetric  relative  difference  between  the 
cell  face  scalar  flux  for  each  cell  width  and  the  corresponding  baseline  solution.  The 
baseline  solution  for  each  combination  was  computed  using  a  cell  width  of  4.88  x 
10-6.  The  observed  rate  of  convergence  for  the  different  spatial  quadratures  for  both 
the  S§  and  DP4  angular  quadratures  agrees  well  with  the  analytical  expectations 
from  Larsen  [19]  (Table  5.2).  The  LD  and  LC  spatial  quadratures  are  of  particular 
interest  because  of  their  high  rate  of  convergence.  Larsen  demonstrated  that  the  global 
discretization  error  for  LD  and  LC  was  0( Ax2)  and  0(Ax3),  respectively.  He  further 
demonstrated  that  LD  and  LC  will  exhibit  a  “superconvergence”  in  the  cell-averaged 
angular  fluxes.  The  rate  of  convergence  in  the  cell-averaged  angular  fluxes  will  be  one 
higher  than  the  global  discretization  error.  Furthermore,  Larsen  demonstrated  that 
the  superconvergence  of  LD  and  LC  also  applies  to  cell-edge  angular  fluxes.  Figures 
5.2  to  5.5  show,  for  each  of  the  four  spatial  quadratures,  the  maximum  value  of  the 
symmetric  relative  difference  along  with  the  convergence  rate  produced  using  a  linear 
fit.  For  LD  and  LC,  some  values  were  omitted  when  the  cell  size  was  small  because 
the  solution  had  converged  to  the  benchmark.  For  optically  thick  cells,  some  values 
were  omitted  because  there  was  additional  sources  of  error,  e.g.  DD  had  negative 
fluxes  for  very  optically  thick  cells. 

When  this  test  was  performed  with  a  double-range  Gauss-Legendre  angular 
quadrature  and  the  total  number  of  ordinates  was  still  8  (DP4),  the  number  of  itera¬ 
tions  required  for  convergence  was  unchanged.  Note  that,  for  optically  thick  cells,  the 
omitted  values  for  the  single-range  angular  quadrature  results  were  closer  to  the  con¬ 
vergence  rate  line  than  the  double-range  angular  quadrature  results.  The  converse  is 
true  for  optically  thin  cells.  This  effect  is  the  manifestation  of  the  performance  differ¬ 
ence  of  the  angular  quadratures  in  optically  thick  versus  optically  thin  cells  mentioned 
previously. 
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Figure  5.2:  Spatial  mesh  refinement  convergence  rate  in 
slab  geometry  using  Distribution  Iteration 
and  diamond  difference. 


70 


Maximum  SRD  of  the  Scalar  Flux  Maximum  SRD  of  the  Scalar  Flux 


Cell  Size  (Mean  Free  Paths) 
(b)  Double-Range,  SC 


Figure  5.3:  Spatial  mesh  refinement  convergence  rate  in 
slab  geometry  using  Distribution  Iteration 
and  step  characteristic. 
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Figure  5.4:  Spatial  mesh  refinement  convergence  rate  in 
slab  geometry  using  Distribution  Iteration 
and  linear  discontinuous. 
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Figure  5.5:  Spatial  mesh  refinement  convergence  rate  in 
slab  geometry  using  Distribution  Iteration 
and  linear  characteristic. 
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Table  5.2:  Convergence  rates  for  the  constant  source, 
single  region  problem  using  the  Distribution 
Iteration  algorithm. 


Convergence  Rate 

Spatial  Quadrature 

Single- Range 

Double-Range 

Analytic 

DD 

2.0 

2.0 

2 

SC 

2.0 

2.0 

2 

LD 

3.0 

2.9 

3 

LC 

3.8 

3.9 

4 

5. 1.1. 4  Tests  5  and  6.  The  two  problems  with  the  left  isotropic  surface 
source  (5  and  6)  were  developed  to  exploit  Ganapol’s  TIEL  benchmark.  Both  the  4 
mfp  and  the  infinite  medium  tests  were  also  used  to  test  left/right  symmetry.  Using  a 
spatial  mesh  of  2048  uniformly  spaced  cells  and  the  S32  and  DPi6  angular  quadratures, 
the  infinite  medium  test  was  compared  with  the  TIEL  benchmark  and  the  scalar 
flux  was  within  the  convergence  tolerance  except  when  too  close  to  the  source.  The 
agreement  of  the  scalar  flux  in  proximity  to  the  isotropic  surface  source  was  dependent 
upon  the  angular  quadrature  order;  the  higher  the  order  the  closer  one  could  approach 
the  source. 


5. 1.1.5  Test  7.  Ideally,  the  isotropic  test  would  also  be  used  with 
the  source  iteration  code,  however  this  was  not  possible  because  PARTISN  does  not 
permit  a  symmetry  boundary  and  an  incident  current  to  share  a  surface.  Instead  I 
developed  a  separate  test  using  a  Lambertian  illumination  of  the  left  boundary  and 
left  and  right  vacuum  boundaries  for  comparing  DI  with  the  source  iteration  code. 
The  change  from  isotropic  surface  source  to  the  Lambertian  illumination  was  due  to 
the  need  to  test  the  implementation  of  that  source  type  within  the  DI  test  code.  Using 
a  spatial  mesh  of  512  uniformly  spaced  cells  and  the  Ss  and  DP4  angular  quadratures 
both  DI  and  the  source  iteration  solutions  agreed  to  within  the  agreement  tolerance. 
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5. 1.1. 6  Tests  8  and  9.  The  periodic  interface  problems  (8  and  9)  are 
used  to  assess  performance  in  regimes  where  there  is  a  repeated  change  in  material 
properties.  These  tests  were  inspired  by  the  periodic  horizontal  interface  problem 
presented  by  Azmy  [5].  Even  though  Azrny’s  test  was  designed  for  XY  geometry,  it 
does  pose  an  interesting  test  case  in  slab  geometry  because  it  provides  a  mean  to  study 
the  how  alternating  material  interfaces  and  deep  penetration  affect  the  performance 
of  DI. 

Test  8  creates  a  strong  discontinuity  in  total  cross-section  between  pure-scattering 
regions.  In  the  high  cross-section  regions  there  is  a  strong  scattering  source  causing  an 
isotropic  angular  flux  being  transported  to  the  adjacent  cells.  The  low  cross-section 
regions  there  is  a  weak  scattering  source,  which  results  in  a  near  streaming  transport 
of  the  inward  fluxes  to  the  adjacent  cells. 

These  tests  were  conducted  by  growing  the  problem  by  adding  a  region  pair , 
which  is  just  the  combination  of  the  two  regions  A  and  B.  Using  a  spatial  mesh  of  64 
uniformly  spaced  cells  for  each  region  and  the  Ss  and  DP4  angular  quadratures,  the 
maximum  value  of  the  symmetric  relative  difference  was  within  the  convergence  tol¬ 
erance  when  using  diffusion  synthetic  acceleration  (DSA);  however,  for  unaccelerated 
source  iteration  and  transport  synthetic  acceleration  (TSA),  the  difference  was  larger 
than  the  convergence  tolerance  (Figure  5.6).  Note  that  PARTISN  turned  off  DSA 
after  four  region  pairs  because  it  was  no  longer  effective.  The  reason  for  the  lack  of 
agreement  between  DI  and  the  solutions  produced  by  SI  is  due  to  false  convergence 
on  the  part  of  SI.  When  DSA  was  enabled,  the  solutions  produced  by  DI  and  SI  were 
in  agreement.  When  DSA  was  disabled  at  five  region  pairs,  the  two  solutions  were  no 
longer  in  agreement. 

Test  9  causes  an  alternating  change  in  the  angular  distribution  of  the  flux. 
The  c  =  1  regions  produce  an  outgoing  isotropic  angular  flux  distribution,  while  the 
c  =  0  regions  produce  an  angular  flux  distribution  that  is  biased  towards  the  X- 
axis.  The  combination  of  the  scattering  regions,  which  effectively  divides  the  flux  of 
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neutrons  into  leftward  and  rightward  flows,  and  the  absorbing  regions  results  in  a 
diminishing  neutron  flux.  This  provides  an  opportunity  to  evaluate  the  performance 
in  deep  penetration  problems.  Using  a  spatial  mesh  of  16  uniformly  spaced  cells  for 
each  region  and  the  Sg  and  DP4  angular  quadratures,  the  agreement  between  DI  and 
PARTISN  was  within  the  convergence  tolerance  to  about  28  (unaccelerated)  to  33 
(DSA)  region  pairs.  Again,  the  lack  of  agreement  between  the  DI  solutions  and  SI 
solutions  is  due  to  false  convergence  by  SI. 


5. 1.1. 7  Test  10.  The  two  region  scattering  ratio  test  provides  insight 
into  the  performance  of  a  transport  code  for  a  variety  of  material  type  interfaces 
and  was  performed  with  both  DI  and  PARTISN.  Furthermore,  the  source  iteration 
code  was  run  with  no  source  acceleration,  with  DSA,  and  with  TSA.  The  preferred 
angular  quadrature  for  an  incident  current  problem  is  double-range  Gauss-Legendre 
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(DPn)  because  of  its  superior  performance  in  solving  for  cell  face  angular  fluxes2  and 
was  used  in  all  cases  except  with  TSA.  PART1SN  does  not  support  the  use  of  DPn 
with  TSA;  therefore  only  the  traditional  S„,  quadrature  was  used  with  TSA.  Using  a 
spatial  mesh  of  16  uniformly  spaced  cells  for  each  region  and  the  Sg  and  DP4  angular 
quadratures,  the  scalar  fluxes  calculated  by  D1  in  all  cases  were  within  the  convergence 
tolerance  of  the  solutions  generated  by  PART1SN. 

5.1.2  X Y  Geometry.  Prins  [28]  developed  and  examined  the  performance 
of  D1  in  XY  geometry.  1  implemented  a  series  of  test  cases  similar  to  the  ones  used 
in  slab  geometry.  The  set  of  test  problems  is  presented  in  the  following  list: 

1.  Pure  Absorber-Single  region  (cr  =  1,  c  =  0,  ge  =  4.7,  32  mfp  by  32  mfp)  and 
vacuum  boundaries; 

2.  Pure  Scatterer-Single  region  (a  =  1,  c  =  1,  qe  =  4.7,  32  mfp  by  32  mfp)  and 
vacuum  boundaries; 

3.  Infinite  Medium-Single  region  (a  =  l,c  =  0.5,  qe  =  4.7,  32  mfp  by  32  mfp), 
symmetry  boundaries; 

4.  Constant  Source-Single  region  (a  =  0.1,  c  =  0.5,  qe  =  4.7,  3.2  mfp  by  32  mfp  by 
32  mfp),  vacuum  boundaries; 

5.  Left  Isotropic  Surface  Source  ( Q  =  0.5),  4  mfp-Single  region  (cr  =  1,  c  =  0.9,  qe  = 
0,  4  mfp  by  4  mfp),  left /top/bottom  symmetry  boundaries  and  right  vacuum 
boundary; 

6.  Left  Isotropic  Surface  Source  ( Q  =  0.5),  32  mfp-Single  region  (a  =  l,c  = 
0.9,  qe  =  0,  32  mfp  by  4  mfp),  left/top/bottom  symmetry  boundaries  and  right 
vacuum  boundary; 

7.  Left  Lambertian  Illumination  ( Q  =  4.7),  4  mfp-Single  Region  (a  —  l,c  — 
0.9,  qe  =  0,  4  mfp  by  4  mfp),  vacuum  boundaries; 

2The  S„  angular  quadrature  can  be  off  by  1%  or  more  for  thin  cells 
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The  only  spatial  quadrature  in  XY  geometry  in  common  between  PARTISN  and  my 
DI  code  was  diamond  difference.  Also,  all  comparisons  were  done  using  the  S4  angular 
quadrature.  The  grid  size  in  all  test  cases  was  set  to  20  by  20;  the  only  exception 
being  the  left  isotropic  surface  source  test  case  which  used  a  40  by  40  for  the  4  mfp 
case  and  a  320  by  40  for  the  32  mfp  case. 

The  use  of  the  diamond  difference  spatial  quadrature  is  problematic  in  XY 
geometry  as  it  is  possible  to  generate  negative  fluxes.  While  there  are  flux  fix  up 
methodologies,  there  are  drawbacks  to  using  them.  1  deliberately  decided  to  use 
diamond  difference  in  the  verification  effort  because  of  the  negative  flux  problem. 
1  wanted  to  evaluate  the  performance  of  DI  when  negative  fluxes  are  generated  as 
this  should  make  it  difficult  for  DI  to  converge.  In  addition  to  diamond  difference,  I 
implemented  the  step  characteristic  and  linear  characteristic  spatial  quadratures. 

I  did  not  perform  any  spatial  mesh  or  angular  quadrature  refinement  testing  in 
the  XY  geometry.  The  structure  of  my  DI  code  in  the  XY  and  XYZ  geometries  was 
designed  to  facilitate  debugging  and  to  provide  a  detailed  look  at  how  the  solution  is 
progressing.  As  a  result,  it  is  neither  memory  nor  computationally  efficient.  Also,  the 
implementation  of  the  linear  discontinuous  spatial  quadrature  is  necessary  to  provide 
more  reliable  results  that  can  be  compared  with  PARTISN. 

The  tests  in  which  both  DI  and  PARTISN  (unaccelerated  and  DSA)  agree  to 
within  the  convergence  tolerance  are: 

•  Pure  Absorber  (Test  1) 

•  Infinite  Medium  (Test  3) 

•  Constant  Source,  Single-Region  (Test  4) 

•  Lambertian  Source,  4mfp,  Vacuum  (Test  7) 

In  the  pure  scatterer  test  (Test  2),  only  DSA  agreed  to  within  the  convergence  toler¬ 
ance.  Unaccelerated  source  iteration  required  586  iterations  and  did  not  agree  with 
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either  DI  or  with  DSA.  This  difference  is  due  to  false  convergence  by  unaccelerated 
SI. 

The  isotropic  surface  source  tests  were  used  to  compare  with  the  TIEL  bench¬ 
mark.  The  matrix  albedo  was  not  implemented  in  XY  nor  XYZ  geometry,  thus  a  32 
mfp  problem  was  used  to  emulate  an  infinite  medium.  When  using  diamond  differ¬ 
ence,  DI  was  unable  to  generate  a  reasonable  solution  due  to  the  negative  flux  effect. 
The  step  characteristic  spatial  quadrature  was  much  closer-the  maximum  symmetric 
relative  difference  was  3.0  x  10-2  for  optical  depths  in  excess  of  1  mfp.  Given  the 
performance  constraints  of  my  DI  code  due  to  the  way  it  was  implemented,  I  did  not 
refine  the  spatial  and  angular  mesh  to  refine  the  solution. 

The  two  region,  scattering  ratio  sweep  was  implemented  in  XY  geometry  as  two 
slabs  with  dimensions  2  mfp  by  4  mfp  (Figure  5.7),  with  vacuum  boundaries  on  all 
four  sides.  A  Lambertian  illumination  was  used  on  the  left  boundary.  Of  the  hundred 
combinations  of  scattering  ratios,  there  were  eight  cases  where  DI  and  DSA  did  not 
agree;  in  three  of  the  eight  cases,  DI  had  failed  to  converge  within  500  iterations.  This 
will  be  discussed  further  in  Section  6.1.2. 


Figure  5.7:  Layout  for  the  two  region  test. 

5.1.3  XYZ  Geometry.  Because  XYZ  geometry  was  an  extension  of  the  XY 
implementation,  a  limited  set  of  test  cases  were  required.  The  focus  of  the  testing 
was  to  ensure  that  Z  axis  currents  and  the  spatial  quadratures  were  implemented 
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correctly.  The  set  of  test  problems  performed  in  XYZ  geometry  is  presented  in  the 
following  list: 

1.  Pure  Absorber-Single  region  (a  =  1,  c  =  0,  qe  =  4.7,  32  mfp  by  32  mfp  by  32 
mfp)  and  vacuum  boundaries; 

2.  Pure  Scatterer-Single  region  (a  —  1,  c  —  1,  qe  —  4.7,  32  mfp  by  32  mfp  by  32 
mfp)  and  vacuum  boundaries; 

3.  ID,  Two  Region,  Scattering  ratio  test  (cr  =  1.0,  qe  =  0,  2  mfp  by  4  mfp  by  4 
mfp),  vacuum  boundaries,  left  incident  Lambertian  illumination  ( Q  =  4.7). 

The  spatial  mesh  in  all  cases  was  a  20  by  20  by  20  uniform  cell  spacing  using  the  S4 
angular  quadrature. 

While  using  DD  with  DI  in  XY  geometry  caused  some  problems  in  converging, 
in  XYZ  geometry  DI  exhibited  exceptionally  poor  performance  when  using  diamond 
difference.  For  problems  with  strong  emission  sources,  such  as  the  pure  absorber 
and  pure  scatterer  tests,  DI  and  PARTISN  converged  to  the  same  solutions.  Using  a 
positive  method,  such  as  step,  DI  is  able  to  converge  for  all  the  test  cases.  My  negative 
flux  fix-up  scheme  was  able  to  fix  many  of  the  cases  were  DI  did  not  converge,  however, 
the  scheme  failed  when  the  scattering  ratio  was  larger  than  0.8. 

5.2  Multigroup  Verification  Tests 

Three  multigroup  test  cases  were  implemented  to  verify  that  the  conventional 
multigroup  method  is  compatible  with  the  DI  method.  The  first  test  was  an  ab¬ 
sorptive  downscatter-only  problem  utilizing  three  energy  groups.  The  group  cross 
sections  are  shown  in  Table  5.3,  where  oa  is  the  absorption  cross-section,  at  is  the 
total  cross-section,  and  <rSg^2  is  the  scattering  cross-section  from  group  g  to  group 
2.  The  downscatter  only  test  is  the  simplest  multigroup  problem  and  it  is  a  useful 
test  for  debugging  purposes.  The  next  test  was  an  absorptive  upscatter  test.  The 
group  cross  sections  are  similar  to  those  shown  in  Table  5.3,  the  only  change  being 
that  the  last  group  has  an  upscatter  contribution  to  the  first  group.  The  third  test 
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was  a  downscatter  only  test  where  the  second  energy  group  was  a  strong  scatterer 
(Table  5.4).  This  test  was  designed  to  have  a  different  inward  angular  partial  current 
distribution  between  adjacent  energy  groups.  The  alternating  pattern  would  test  how 
well  DI  is  able  to  adapt  to  the  change  in  the  inward  angular  current  distribution. 

In  slab  geometry,  all  three  test  cases  consisted  of  a  32  mfp  uniform  region  with 
vacuum  boundaries  and  a  256  uniform  spatial  cell  mesh.  Both  the  Ss  and  DP4  angular 
quadratures  were  used.  In  XY  geometry  a  32  mfp  by  32  mfp  uniform  region  and  20 
by  20  uniform  spatial  cell  mesh  was  used.  The  S4  angular  quadrature  was  used  for 
XY  geometry.  For  the  upscatter  test,  the  convergence  test  was  taking  the  maximum 
symmetric  relative  difference  between  the  current  and  previous  iteration  for  each 
energy  group.  The  convergence  tolerance  was  set  to  10-7. 

DI  and  PARTISN  (unaccelerated  and  TSA)  both  agreed  using  all  combinations 
of  angular  quadratures  and  supported  spatial  quadratures  (LD  and  DD  in  slab  geom¬ 
etry  and  DD  in  XY  geometry).  PARTISN  did  not  agree  to  within  the  convergence 
tolerance  on  the  “Three  Group  Absorber  Upscatter  Test”  when  using  DSA  in  both 
slab  and  XY  geometries.  When  using  DSA,  the  maximum  value  of  the  symmetric 
relative  difference  for  the  different  combinations  of  spatial  and  angular  quadratures 
was  7  x  10”6  to  9  x  1CU6,  just  beyond  the  agreement  requirement  of  10-6. 

5.3  Time  Dependent  Verification  Tests 

Three  time-dependent  external  sources,  qe,  were  used  for  testing.  The  first 
source  was  a  uniform  isotropic  source  that  is  initially  off  and  instantaneously  turns 


Table  5.3:  Group  cross  sections  use  for  the  downscatter 
test. 


Energy  Group 

cra 

crt 

crss_i 

aSg^2 

aSg^3 

1 

0.9 

1.0 

0.0 

0.1 

0.0 

2 

0.9 

1.0 

0.0 

0.0 

0.1 

3 

1.0 

1.0 

0.0 

0.0 

0.0 
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Table  5.4:  Group  cross  sections  use  for  the  alternating 
test. 


Energy  Group 

(7  a 

i A 

aSg^l 

aSg^2 

°Sg_  3 

1 

0.9 

1.0 

0.0 

0.1 

0.0 

2 

0.0 

1.0 

0.0 

0.9 

0.1 

3 

1.0 

1.0 

0.0 

0.0 

0.0 

on  and  then  stays  constant  in  time  (Step  Source).  The  “Step  Source”  test  is  the  most 
rudimentary  time  dependent  test  and  is  useful  in  debugging.  This  test  was  also  used 
to  verify  that  the  solution  was  approaching  the  steady-state  solution.  The  second 
and  third  sources  were  a  “Ramp  Up  and  Hold”  and  a  “On  and  Ramp  Down.”  The 
combination  of  these  three  tests  would  test  the  performance  of  the  DI  method  for  the 
common  time-dependent  source  types.  In  all  test  case  a  single  energy  group  was  used 
with  a  velocity  of  10_3cm/s. 

Comparing  PARTISN  and  DI  on  a  time  step  basis  was  not  feasible  because  of 
the  adaptive  time  step  algorithm  that  PARTISN  uses.  I  decided  that  implementing 
PARTISN’s  adaptive  time  step  algorithm  was  not  useful  nor  was  modifying  PARTISN 
to  create  a  hie  containing  the  time  steps.  Instead,  I  had  DI  use  a  smaller  time 
step,  10~4  seconds,  and  interpolated  between  time  steps  in  order  to  compare  with 
PARTISN.  The  basis  for  my  decision  was  that  PARTISN’s  utility  as  an  accurate 
benchmark  is  limited  by  its  inability  to  output  the  more  accurate  extrapolated  fluxes 
(see  Appendix  A).  Furthermore,  the  focus  of  my  research  was  to  demonstrate  that  DI 
did  not  have  problems  converging  when  used  with  very  optically  thick  spatial  cells. 
The  time-independent  testing  demonstrated  that  DI  produces  reliable  results  when 
the  spatial  cell  is  the  appropriate  size  for  the  spatial  quadrature. 

The  synthetic  cross  section  for  the  test  cases  performed  by  DI  was  10'  cnr1. 
For  PARTISN,  the  synthetic  cross  section  varied  between  107  cnr1  and  104  cnr1. 
The  constant  source  test  was  also  used  to  test  a  range  of  time  steps,  ranging  from  10 
seconds  to  10”5  seconds. 
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Table  5.5  shows  the  maximum  and  mean  values  for  the  symmetric  relative  dif¬ 
ference  between  the  D1  and  PART1SN  solutions  taken  at  each  PART1SN  time  step. 
The  interpolation  was  performed  using  Mathematica’s  interpolation  function,  which 
performs  a  piecewise  polynomial  fit  of  order  3.  This  test  was  performed  when  D1  was 
set  to  output  the  midpoint  scalar  flux  and  the  more  accurate  end  of  time  step  scalar 
flux.  There  was  no  significant  difference  when  the  time  step  was  varied. 

For  XY  geometry,  a  corner  source  problem  was  used.  The  problem  was  6  cm  by 
6  cm  with  an  uniform  isotropic  source  region  in  the  lower  left  corner  with  dimensions 
1  cm  by  1  cm.  The  source  region  material  had  a  cross  section  of  0.1cm-1  and  a 
scattering  ratio  of  0.5.  The  rest  of  the  problem  was  a  uniform  material  with  a  cross 
section  of  1.0cm-1  and  a  scattering  ratio  of  0.5.  A  20  by  20  spatial  mesh  was  used 
and  the  angular  quadrature  was  S4. 

In  XY  geometry,  comparing  the  results  generated  by  PARTISN  and  DI  was 
problematic.  When  using  DD  as  the  spatial  quadrature  a  large  number  of  negative 
scalar  fluxes  will  occur  outside  the  source  region.  PARTISN  forces  a  negative  flux  fix¬ 
up  for  time  dependent  problems,  thus  it  is  not  possible  to  compare  my  DI  code  and 
PARTISN  without  any  fix-up.  Because  I  did  not  implement  the  same  negative  flux  fix¬ 
up,  there  is  a  substantial  difference,  a  symmetric  relative  difference  of  approximately 
0.1,  in  the  results  generated  by  DI  and  PARTISN. 


Table  5.5:  Maximum  and  mean  symmetric  relative  dif¬ 
ference  between  DI  and  PARTISN  for  the 
time  dependent  tests  in  slab  geometry. 


Test 

Maximum 

Mean 

Constant  Source 

1.0  x  10-4 

2.2  x  10-6 

On  and  Ramp  Down 

6.4  x  10-4 

6.3  x  10-4 

Ramp  Up  and  Hold 

6.0  x  10-4 

8.4  x  10-5 
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VI.  Performance  Analysis 


The  verification  testing  that  I  performed  demonstrated  that  DI  converges  to  the 
benchmark  solution  in  most  of  the  test  cases.  An  assessment  of  how  the  Df  algorithm 
performed  in  reaching  the  solution  is  necessary  in  order  to  judge  if  Dl  is  a  suitable 
replacement  for  source  iteration.  The  key  performance  metric  in  this  assessment  is 
the  number  of  iterations  required  for  convergence.  Because  of  the  different  perfor¬ 
mance  trade-offs  made  in  my  Dl  code,  which  is  designed  to  support  research,  versus 
PART1SN,  which  is  designed  to  be  a  production  code,  a  run  time  or  processor  time 
comparison  is  not  a  reliable  measure  of  performance. 

In  the  case  of  PARTISN,  the  number  of  inner  iterations  is  the  measure  of  its 
performance.  For  DI,  the  number  of  calls  to  the  partial  current  solver  is  the  measure 
of  its  performance.  In  addition,  an  examination  of  when  DI  did  not  converge,  or 
converged  to  a  different  solution,  is  part  of  the  qualitative  assessment  of  the  overall 
performance  of  the  DI  method. 

6.1  Time  Independent  Problems 

All  the  test  cases  referenced  in  this  chapter  are  described  in  Chapter  V. 

6.1.1  Slab  Geometry.  A  summary  of  the  number  of  iterations  required 
for  convergence  is  shown  in  Table  6.1.  In  all  cases  DI  was  able  to  converge  to  the 
benchmark  solution  in  the  same  number  or  fewer  iterations  as  PARTISN  using  DSA. 
Furthermore,  with  few  exceptions,  the  agreement  between  the  scalar  fluxes  produced 
by  DI  and  PARTISN  was  within  the  convergence  tolerance.  The  few  exceptions  where 
the  two  methods  did  not  agree  are  presented  later  in  the  text.  The  isotropic  source 
tests  (5  and  6)  was  not  performed  with  PARTISN  because  the  source  could  not  be 
accurately  represented.  The  ranges  given  for  tests  8  and  9  reflect  the  number  of 
iterations  from  the  smallest  problem  to  the  largest  problem. 

6. 1.1.1  Tests  1  and  2.  The  single  region  pure  absorber  test  (1)  demon¬ 
strates  that  there  is  no  performance  advantage  between  DI  and  source  iteration  for 
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Table  6.1:  Number  of  iterations  required  for  convergence 
in  the  time  independent, slab  geometry  tests. 


Benchmark 

DI 

Source  Iteration 

None  DSA 

TSA 

1.  Pure  Absorber 

2 

2 

2 

2 

2.  Pure  Scatterer 

7 

103 

10 

25 

3.  Infinite  Medium 

1 

27 

10 

N/A 

4.  Constant  Source,  Single-Region 

6 

20 

7 

7 

5.  Isotropic  Source,  4  mfp 

11 

N/A 

N/A 

N/A 

6.  Isotropic  Source,  Infinite  Medium 

11 

N/A 

N/A 

N/A 

7.  Lambertian  Source,  4  mfp,  Vacuum 

8 

75 

15 

18 

8.  Periodic  Interface,  Cross-sections 

9  -  12 

540  -  27891 

22  -  41684  123 

-  6605 

9.  Periodic  Interface,  Scattering  ratio 

8  -  17 

33  -  98 

11  -  22 

12  -  29 

strong  absorbing  problems.  As  the  scattering  ratio  increases  (Test  4  and  becomes 
a  pnre  scattering  problem  (Test  2),  the  performance  of  unaccelerated  source  itera¬ 
tion  drops.  There  is  no  performance  advantage  between  D1  and  accelerated  source 
iteration,  either  DSA  or  TSA,  for  problems  that  feature  uniform  material  properties. 

6. 1.1. 2  Test  3.  This  test  problem  is  solved  in  one  iteration  by  D1 
because  the  initial  inward  current  distribution  is  isotropic,  which  is  correct  distribu¬ 
tion  for  this  problem.  If  the  starting  inward  current  distribution  is  set  randomly,  D1 
requires  eight  iterations  to  converge.  This  test  suggests  a  useful  performance  enhance¬ 
ment  for  D1  by  using  an  approximation  for  the  inward  current  distribution.  This  is 
particularly  useful  in  a  time-dependent  code  because  in  a  slowly-varying  problem  the 
current  distribution  will  not  change  rapidly  between  time  steps.  D1  performed  slightly 
better  than  DSA  with  this  problem. 

6. 1.1.3  Test  4-  Generally,  all  four  spatial  quadratures  converged 
within  6  iterations,  the  only  exceptions  being  at  the  coarse  grid  spacings  where  LC, 
SC,  and  DD  required  5  iterations  (SC  required  only  4  iterations  at  the  coarsest  mesh 
of  8  cells).  As  previously  noted,  there  was  no  performance  advantage  between  D1  and 
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accelerated  source  iteration.  DI  performed  much  better  than  unaccelerated  source 
iteration. 


6. 1.1. 4  Test  7.  The  Lambertian  illumination  problem  exhibited  a  clear 
performance  advantage  between  Dl  and  accelerated  source  iteration.  DI  was  able  to 
converge  in  half  as  many  iterations  as  DSA.  Figures  6.1  through  6.4  show  the  number 
of  iterations  required  for  convergence  as  the  scattering  varied  from  0  to  1  when  using 
DSA,  TSA,  unaccelerated  PARTISN  and  DI.  Two  test  scenarios  are  used:  One  is  a 
uniform  source  and  the  other  is  a  Lambertian  illumination.  For  DSA,  Figure  6.1, 
there  is  a  difference  in  the  number  of  iterations  between  the  two  tests  when  c  >  0.2. 
For  TSA,  Figure  6.2,  the  difference  in  the  number  of  iterations  does  not  occur  until 
c  >  0.7.  Unaccelerated  source  iteration  (Figure  6.3)  is  similar  to  TSA  in  behavior. 
While  DI  (Figure  6.4)  does  exhibit  the  difference  in  the  number  of  iterations,  there  is 
no  growth  as  c  — >  1. 
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Figure  6.1:  Number  of  iterations  needed  for  Conver¬ 
gence  using  PARTISN  with  DSA,  LD  and 
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Figure  6.2:  Number  of  iterations  needed  for  conver¬ 
gence  using  PARTISN  with  TSA,  LD  and 
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Figure  6.3:  Number  of  Iterations  Needed  for  Conver¬ 
gence  Using  Unaccelerated  PARTISN,  LD 
and  S8. 
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The  Lambertian  illumination  test  is  more  difficult  to  solve  than  a  uniform  source 
because  of  coupling  of  the  flux  between  the  left  and  right  halves  of  the  problem  and 
the  asymmetry  in  the  flux.  When  the  angular  flux  is  perturbed  in  the  left  hand  side, 
it  causes  a  perturbation  in  the  right  hand  side,  which  then  couples  back  to  the  left 
hand  side  (provided  c  >  0).  The  presence  of  an  embedded  neutron  source  masks  the 
effect  of  perturbation  and,  hence,  allows  the  problem  to  converge  faster. 

6. 1.1.5  Tests  8  and  9.  For  Test  8,  both  unaccelerated  source  iteration 
and  TSA  required  more  iterations  than  DI  to  converge.  DSA  had  good  performance, 
requiring  only  22  to  27  iterations  when  the  problem  was  smaller  than  5  region  pairs. 
When  the  problem  was  5  region  pairs  and  larger,  DSA  was  not  providing  any  benefit 
and  was  turned  off  by  PARTISN,  thus  effectively  becoming  an  unaccelerated  method. 

DI  and  PARTISN  did  not  agree  to  within  the  convergence  tolerance  when  PAR¬ 
TISN  required  many  iterations  to  converge.  For  example,  TSA  and  DI  agreed  in  the 
one  region  pair  case  (123  iterations)  and  did  not  agree  in  the  two  region  pair  case 
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Figure  6.4:  Number  of  Iterations  Needed  for  Conver¬ 
gence  Using  DI  With  LD  and  Ss- 


(386  iterations).  The  source  of  the  discrepancy  is  difficult  to  fully  ascertain  without 
instrumenting  PARTISN.  The  most  likely  cause  is  a  combination  of  false  convergence. 

For  Test  9,  DI  required  fewer  iterations  to  converge  (Figure  6.5)  than  PAR¬ 
TISN.  As  the  problem  got  larger  than  30  region  pairs,  PARTISN  would  not  perform 
any  additional  iterations.  The  reason  for  the  difference  between  the  solutions  pro¬ 
duced  by  DI  and  PARTISN  is  due  to  false  convergence  on  the  part  of  PARTISN.  The 
convergence  test  that  PARTISN  performs  is 


max 

i,9 


<  e, 


(6.1) 


where  4>\  is  the  scalar  flux  at  iteration  l  at  mesh  point  i.  group  g  and  e  is  the 
convergence  tolerance.  At  30  region  pairs,  the  scalar  flux  at  the  right  edge  of  the 
problem  is  2.08  x  10~24.  Recalling  (1.24)  and  the  problem  of  small  contributions 
as  the  number  of  scatters  gets  large,  PARTISN  is  exhibiting  that  behavior.  With 
such  a  small  flux,  the  contribution  from  neutrons  that  have  over  95  collisions  (using 
unaccelerated  source  iteration  as  the  measure  of  scattering  contributions)  to  the  scalar 
flux  is  smaller  than  the  convergences  tolerance. 

Overall,  agreement  between  DI  and  PARTISN  for  these  two  tests  is  acceptable 
given  the  limitations  of  the  source  iteration  method.  Further  research  is  needed  to 
verify  the  solutions  generated  by  DI  beyond  the  regime  where  PARTISN  generates 
acceptable  results. 


6. 1.1. 6  Test  10.  Tables  6.2  through  6.5  show  the  number  of  iterations 
required  for  convergence  for  test  10  (the  Two  Region  Scattering  Ratio  Sweep  Test).  I 
also  ran  the  test  using  the  LD  spatial  quadrature  and  observed  no  significant  difference 
in  the  number  of  iterations.  The  number  of  iterations  required  for  convergence  was 
essentially  the  same  for  both  the  DPn  and  Sn  angular  quadratures.  In  all  cases,  DI 
converged  in  fewer  iterations,  particularly  as  c  — »  1. 
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Figure  6.5:  Number  of  Iterations  Needed  for  Conver¬ 
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The  results  shown  in  tables  6.3  and  6.5  are  symmetric  relative  to  the  pnre 
absorber  to  the  pure  scatter  diagonal  (a  difference  of  one  iteration  is  ignored).  Table 
6.2  has  a  slight  asymmetry  when  the  absorbing  region  is  on  the  left  versus  the  right. 
When  the  absorbing  region  is  on  the  left  and  the  scattering  region  is  on  the  right, 
there  are  two  isotropic  (or  nearly  isotropic)  flows  into  the  left  region.  This  forces 
DI  to  transition  from  isotropic  to  an  anisotropic  flow  from  two  different  directions. 
Conversely,  when  the  absorbing  region  is  on  the  right  and  the  scattering  region  is  on 
the  left,  there  is  only  one  isotropic  flow  into  the  right  region. 

The  DSA  case  (Table  6.4)  also  has  an  asymmetry  relative  to  the  diagonal.  When 
the  scattering  region  is  on  the  left  and  the  absorber  is  on  the  right,  the  problem  is 
similar  to  the  Left  Lambertian  Illumination  test  (Test  7).  When  the  right  region  is 
a  pnre  absorber,  the  problem  is  equivalent  to  a  2  mfp  Left  Lambertian  Illumination 
test  (which  converges  in  the  same  number  of  iterations).  As  the  scattering  ratio  of  the 
right  region  increases,  the  problem  evolves  into  a  4  mfp  Left  Lambertian  Illumination 
test.  When  the  absorbing  region  is  on  the  left  and  the  scattering  region  is  on  the  right, 
the  left  region  is  quickly  solved  and  the  right  region  has  an  anisotropic  illumination 
on  the  left. 

6.1.2  X Y  Geometry.  A  summary  of  the  number  of  iterations  required  for 
convergence  is  shown  in  Table  6.6.  The  performance  that  DI  exhibited  when  using 
the  diamond-difference  spatial  quadrature  did  not  reveal  any  significant  performance 
advantage;  in  fact,  in  some  cases  DI  required  more  iterations.  For  example,  in  the 
case  of  the  pnre  scatterer  test,  DI  needed  more  iterations  to  converge  than  DSA  when 
using  DD;  however,  with  SC  and  LC,  DI  required  only  9  iterations  to  converge. 

The  infinite  medium  test,  as  in  the  slab  geometry  case,  was  able  to  converge 
within  one  iteration  when  using  an  isotopic  distribution  as  the  initial  condition. 
Changing  DI  to  use  a  random  initialization  of  the  inward  current  distribution  pre¬ 
vents  convergence  when  using  DD.  If  SC  is  used  instead  of  DD,  DI  is  able  to  converge 
within  17  iterations. 
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Table  6.2:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  DI,  diamond  differ¬ 
ence,  DP4  double-range  Gauss-Legendre  an¬ 
gular  quadrature. 


0.0 

0.1 

Left  Region  Scattering  Ratio 

0.2  0.3  0.4  0.5  0.6  0.7  0.8 

0.9 

1.0 

.2 

0.0 

2 

6 

6 

7 

7 

7 

7 

8 

8 

8 

7 

c5 

0.1 

6 

6 

6 

7 

7 

7 

8 

8 

8 

8 

8 

bO 

£3 

0.2 

7 

7 

7 

7 

7 

7 

8 

8 

8 

8 

8 

'C 

CP 

0.3 

7 

7 

7 

7 

7 

7 

8 

8 

8 

8 

8 

0.4 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

O 

co 

0.5 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

£ 

0 

0.6 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

8 

‘3d 

0D 

0.7 

9 

9 

9 

8 

8 

8 

8 

8 

8 

8 

8 

0.8 

9 

9 

9 

9 

9 

9 

9 

8 

8 

8 

8 

4=1 

bC 

0.9 

9 

9 

9 

9 

9 

9 

9 

8 

8 

8 

8 

s 

1.0 

9 

9 

9 

9 

9 

9 

9 

8 

8 

8 

7 

Table  6.3:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  SI,  no  source  acceler¬ 
ation,  diamond  difference,  DP4  double-range 
Gauss-Legendre  angular  quadrature. 


Left  Region  Scattering  Ratio 


0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

.0 

0.0 

2 

8 

11 

13 

16 

19 

23 

28 

35 

46 

63 

cb 

0.1 

8 

9 

11 

13 

16 

19 

23 

28 

36 

46 

64 

bO 

0.2 

11 

11 

12 

14 

16 

20 

24 

29 

36 

47 

65 

'C 

CD 

0.3 

13 

13 

14 

15 

17 

20 

24 

29 

36 

47 

66 

+=> 

0.4 

16 

16 

17 

17 

19 

21 

25 

30 

37 

48 

67 

O 

CO 

0.5 

19 

19 

20 

20 

21 

23 

27 

32 

39 

50 

70 

£ 

O 

0.6 

23 

23 

24 

24 

25 

27 

29 

34 

41 

52 

73 

‘3d 

QD 

0.7 

28 

28 

29 

29 

30 

32 

34 

38 

45 

56 

78 

0.8 

35 

35 

36 

37 

38 

39 

41 

45 

51 

63 

87 

=1 

_bC 

0.9 

45 

46 

47 

47 

49 

50 

53 

56 

63 

75 

103 

2 

1.0 

62 

63 

64 

66 

68 

70 

73 

78 

87 

103 

141 

92 


Table  6.4:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  SI,  DSA,  diamond 
difference,  DP4  double-range  Gauss-Legendre 
angular  quadrature. 


Left  Region  Scattering  Ratio 


0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

.2 

0.0 

2 

6 

7 

8 

8 

9 

10 

11 

12 

13 

14 

c5 

0.1 

5 

6 

7 

8 

8 

9 

10 

11 

12 

13 

14 

bO 

£3 

0.2 

6 

6 

7 

7 

8 

9 

10 

11 

12 

13 

14 

'C 

CD 

0.3 

7 

7 

7 

7 

8 

9 

10 

11 

12 

13 

15 

4^> 

0.4 

7 

7 

7 

7 

8 

9 

10 

11 

12 

13 

15 

O 

m 

0.5 

8 

8 

8 

8 

8 

9 

10 

11 

12 

14 

15 

£ 

0 

0.6 

8 

8 

8 

8 

8 

9 

10 

11 

12 

14 

15 

‘3d 

0D 

0.7 

9 

9 

9 

8 

9 

9 

10 

11 

12 

14 

16 

0.8 

10 

9 

9 

9 

9 

9 

10 

11 

13 

14 

16 

4-^ 

4=1 

bC 

0.9 

10 

10 

10 

10 

9 

9 

10 

11 

13 

15 

17 

s 

1.0 

11 

11 

11 

11 

11 

10 

10 

11 

13 

15 

18 

Table  6.5:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  SI,  TSA,  diamond 
difference,  DP4  double-range  Gauss-Legendre 
angular  quadrature. 


Left  Region 

Scattering  Ratio 

0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

_o 

0.0 

2 

5 

6 

7 

7 

8 

8 

9 

10 

13 

17 

cb 

0.1 

5 

5 

6 

7 

7 

8 

8 

9 

10 

13 

17 

bO 

0.2 

6 

6 

6 

7 

7 

8 

8 

9 

10 

13 

17 

'C 

CD 

0.3 

7 

6 

6 

7 

7 

8 

8 

9 

10 

13 

17 

4^ 

4^ 

0.4 

7 

7 

7 

7 

7 

8 

8 

9 

10 

13 

18 

O 

m 

0.5 

8 

8 

7 

7 

8 

8 

8 

9 

11 

13 

18 

£ 

0 

0.6 

8 

8 

8 

8 

8 

8 

9 

9 

11 

14 

19 

‘3d 

QD 

0.7 

9 

9 

9 

9 

9 

9 

9 

10 

11 

14 

19 

0.8 

10 

10 

10 

10 

10 

11 

11 

11 

13 

15 

21 

4— ’ 

43 

_bC 

0.9 

13 

13 

13 

13 

13 

13 

14 

14 

15 

18 

25 

£ 

1.0 

17 

17 

17 

17 

18 

18 

19 

19 

21 

25 

33 
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In  the  two  region  scattering  ratio  sweep  test,  the  performance  exhibited  by  DI 
when  using  the  diamond-difference  spatial  quadrature  (Table  6.7)  was  significantly 
worse  than  PARTISN  using  DSA  (Table  6.11).  For  some  combinations  of  scattering 
ratio,  DI  was  unable  to  converge  in  fewer  than  500  iterations  (indicated  by  “NC”  in 
the  table).  When  compared  to  unaccelerated  source  iteration  (Table  6.10),  DI  with 
diamond  difference  performs  better  in  most  cases. 

The  eight  cases  where  DI  and  PARTISN  (indicated  in  bold  in  Table  6.7)  do 
not  have  any  apparent  pattern  other  than  they  all  occur  when  the  left  region  has 
a  scattering  ratio  less  than  0.7  and  in  five  of  the  eight  cases  the  left  region  had  a 
scattering  ratio  less  than  0.5.  When  there  are  no  embedded  sources  in  a  cell  and  the 
scattering  ratio  is  low,  negative  fluxes  are  more  likely  to  occur.  The  negative  fluxes 
will  result  in  a  variety  of  unphysical  artifacts,  e.g.  strong  axial  flows  and  oscillations 
in  the  angular  distribution,  which  is  the  most  likely  cause  for  the  differences  in  the 
two  solutions.  Interestingly  enough,  only  three  of  the  four  cases  where  DI  did  not 
converge  within  500  iterations  are  in  this  set  of  eight;  one  of  the  non-convergent  cases 
was  able  to  produce  the  same  solution  as  PARTISN.  In  all  eight  cases,  the  area  of  the 
problem  where  the  two  solutions  failed  to  agree  was  always  along  the  far  right  edge. 
The  angular  flux  in  the  right  half  of  the  problem  is  much  smaller  than  the  angular 
flux  in  the  left  half  because  the  only  source  of  neutrons  is  from  the  illumination  on 
the  left  boundary,  thus  negative  fluxes  are  more  likely  to  occur. 

When  the  negative  flux  fix-up  scheme  was  used,  the  irregular  pattern  of  iteration 
counts  vanished  (Table  6.8).  There  was  an  increase  in  iteration  count  when  c  >  0.8 
relative  to  PARTISN.  A  direct  comparison  of  the  scalar  fluxes  generated  by  DI  and 
PARTISN  is  not  possible  because  the  difference  in  fix-up  schemes  will  cause  them  to 
converge  to  different  solutions.  Instead,  I  compared  the  DI  and  PARTISN  results  to 
the  results  produced  by  DI  using  the  step  characteristic  spatial  quadrature.  Because 
step  characteristic  is  non-negative  and  has  second-order  convergence,  this  comparison 
would  highlight  the  effect  of  the  fix-up  on  the  solution.  The  magnitude  of  the  variation 
in  the  two  diamond  difference  solutions  relative  to  the  step  characteristic  solution  was 
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approximately  the  same.  Thus,  the  solution  produced  by  my  fix-up  scheme  was  no 
worse  than  the  solution  produced  by  PARTISN  with  its  fix-up  scheme. 

When  using  a  positive  spatial  quadrature  such  as  step  characteristic,  the  per¬ 
formance  of  DI  improves  dramatically  (Table  6.9).  This  behavior  further  highlights 
the  problem  of  using  a  non-positive  spatial  quadrature  with  distribution  iteration. 

6.1.3  XYZ  Geometry.  The  two  region  sweep  problem  (Test  10)  proved  to 
be  very  challenging  for  DI.  For  virtually  all  of  the  combinations  of  scattering  ratios, 
DI  was  unable  to  converge  within  50  iterations  when  using  diamond-difference,  even 
after  the  convergence  tolerance  was  changed  to  1.0  x  10”5.  The  flux  fix-up  scheme  was 
able  to  correct  most  of  the  non-convergent  cases,  however,  when  c  >  0.7  the  fix-up 
scheme  lost  effectiveness  in  reducing  the  iteration  count. 

Table  6.12  summarizes  the  number  of  iterations  required  for  convergence.  The 
tests  where  DI  was  unable  to  converge  using  DD  are  noted  with  “NC.”  When  a 
positive  spatial  quadrature  was  used,  DI  was  able  to  converge. 

6.2  Multigroup  Verification  Tests 

Table  6.13  summarizes  the  number  of  iterations  that  each  method  required  to 
reach  the  convergence  tolerance  in  slab  geometry.  DI  performed  as  well,  or  better, 
than  either  unaccelerated  source  iteration  or  TSA.  DSA  converged  more  rapidly  in 
the  “Three  Group  Absorber,  Upscatter”  test;  however,  DSA  did  not  converge  to  the 
same  solution  as  DI,  unaccelerated  source  iteration  and  TSA.  The  PARTISN  log  file 
did  not  give  any  indication  that  it  encountered  any  problems  and  had  to  disable  DSA. 
Without  performing  a  detailed  analysis  of  the  PARTISN  code,  it  is  not  clear  why  DSA 
did  not  converge  to  the  same  solution. 

The  results  for  XY  geometry  was  similar  to  slab  geometry,  including  the  differ¬ 
ence  in  the  “Three  Group  Absorber,  Upscatter”  test.  For  the  downscatter  tests  the 
number  of  iterations  was  the  same  as  in  the  slab  geometry  case.  In  XY  geometry,  my 
DI  code  initializes  the  angular  inward  current  distribution  for  every  outer  iteration, 
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Table  6.6:  Number  of  iterations  required  for  convergence 
in  the  time  independent,  XY  geometry  tests. 


DI 

Source  Iteration  with  DD 

Test  Problem 

DD 

SC 

LInaccelerated 

DSA 

Pure  Absorber 

2 

2 

2 

2 

Pure  Scatterer 

16 

4 

586 

11 

Infinite  Medium 

1 

1 

139 

106 

Constant  Source,  Single-Region 

6 

5 

18 

7 

Isotropic  Source,  4mfp 

31 

23 

N/A 

N/A 

Isotropic  Source,  32mfp 

NC 

65 

N/A 

N/A 

Lambertian  Source,  4mfp,  Vacuum 

18 

8 

22 

9 

Table  6.7:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  DI,  diamond  differ¬ 
ence,  S4  single-range  Gauss-Legendre  angular 
quadrature  in  XY  geometry. 


Left  Region  Scattering  Ratio 


0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

.2 

0.0 

2 

5 

6 

8 

7 

8 

9 

10 

12 

11 

15 

cS 

0.1 

8 

51 

14 

9 

7 

9 

10 

10 

12 

11 

15 

bO 

0.2 

8 

8 

9 

13 

8 

14 

NC 

10 

11 

11 

13 

'C 

CD 

0.3 

12 

13 

8 

8 

78 

11 

12 

13 

12 

13 

14 

4-^ 

0.4 

13 

10 

12 

25 

15 

NC 

11 

10 

11 

13 

16 

O 

0.5 

11 

149 

15 

22 

11 

18 

10 

15 

13 

13 

16 

£ 

0 

0.6 

16 

13 

25 

19 

10 

12 

20 

11 

13 

14 

15 

'So 

CD 

0.7 

56 

15 

17 

12 

14 

12 

18 

11 

24 

14 

16 

Pi 

0.8 

17 

15 

14 

17 

14 

12 

14 

36 

14 

15 

18 

4-2 

4=: 

_  bp 

0.9 

20 

14 

18 

16 

22 

NC 

14 

NC 

14 

17 

18 

s 

1.0 

19 

15 

23 

27 

24 

17 

15 

17 

17 

20 

21 
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Table  6.8:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  DI,  diamond  differ¬ 
ence,  S4  single-range  Gauss-Legendre  angular 
quadrature  in  XY  geometry  and  a  negative 
flux  fix-up. 


0.0 

0.1 

Left  Region  Scattering  Ratio 

0.2  0.3  0.4  0.5  0.6  0.7  0.8 

0.9 

1.0 

.2 

0.0 

2 

5 

7 

6 

7 

7 

9 

10 

11 

11 

13 

aj 

Ptf 

0.1 

5 

5 

6 

7 

8 

8 

9 

9 

11 

11 

13 

hX) 

0.2 

6 

6 

6 

7 

7 

8 

8 

10 

10 

11 

13 

'C 

CD 

0.3 

6 

7 

7 

7 

8 

8 

9 

11 

11 

11 

14 

0.4 

8 

7 

8 

8 

8 

9 

10 

10 

11 

12 

14 

CD 

c n 

0.5 

8 

7 

8 

8 

9 

9 

10 

10 

12 

12 

14 

3 

0 

0.6 

9 

9 

9 

9 

9 

9 

13 

11 

12 

13 

14 

‘5o 

CD 

0.7 

10 

9 

10 

10 

12 

11 

11 

11 

12 

14 

16 

0.8 

10 

11 

11 

11 

11 

12 

13 

13 

13 

15 

16 

2 

bp 

0.9 

12 

11 

12 

13 

13 

13 

13 

15 

15 

15 

16 

s 

1.0 

13 

13 

14 

14 

15 

16 

16 

19 

16 

19 

21 

Table  6.9:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  DI,  step  character¬ 
istic,  S4  single-range  Gauss-Legendre  angular 
quadrature  in  XY  geometry. 


0.0 

0.1 

0.2 

Left  Region  Scattering  Ratio 
:  0.3  0.4  0.5  0.6  0.7  0.8 

0.9 

1.0 

.2 

0.0 

2 

4 

5 

5 

6 

6 

7 

7 

8 

8 

9 

cd 

0.1 

5 

5 

5 

6 

6 

6 

7 

7 

8 

9 

9 

bO 

0.2 

5 

5 

6 

6 

6 

7 

7 

8 

8 

9 

9 

Th 

CD 

0.3 

6 

6 

6 

6 

6 

7 

7 

8 

8 

9 

10 

0.4 

6 

6 

6 

7 

7 

7 

8 

8 

8 

9 

10 

CD 

CO 

0.5 

7 

7 

7 

7 

7 

8 

8 

8 

9 

9 

10 

£ 

O 

0.6 

7 

7 

7 

8 

8 

8 

8 

8 

9 

10 

10 

‘5o 

CD 

0.7 

8 

8 

8 

8 

8 

8 

9 

9 

10 

10 

10 

0.8 

9 

9 

9 

9 

9 

9 

9 

10 

10 

10 

11 

2 

bp 

0.9 

9 

9 

10 

10 

10 

10 

10 

10 

10 

11 

12 

2 

1.0 

10 

10 

10 

10 

11 

11 

11 

11 

12 

12 

13 
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Table  6.10:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  SI,  no  acceleration, 
diamond  difference,  S4  single-range  Gauss- 
Legendre  angular  quadrature  in  XY  geome¬ 
try. 


Left  Region  Scattering  Ratio 


0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

.2 

0.0 

2 

9 

11 

15 

16 

19 

22 

26 

32 

39 

51 

aj 

Ptf 

0.1 

9 

10 

12 

14 

16 

19 

22 

26 

32 

39 

51 

hX) 

0.2 

12 

12 

13 

14 

16 

19 

22 

27 

32 

40 

51 

'C 

CD 

0.3 

14 

14 

14 

15 

17 

20 

23 

27 

33 

40 

52 

0.4 

16 

16 

16 

17 

18 

20 

23 

28 

33 

41 

53 

CD 

c n 

0.5 

19 

19 

19 

20 

20 

22 

25 

29 

34 

42 

54 

£ 

0 

0.6 

22 

22 

22 

23 

24 

25 

27 

30 

35 

43 

56 

‘5o 

CD 

0.7 

26 

26 

27 

27 

28 

29 

30 

33 

38 

46 

58 

0.8 

32 

32 

32 

33 

33 

34 

36 

38 

42 

49 

62 

'bp 

0.9 

39 

40 

40 

40 

41 

42 

43 

46 

49 

56 

69 

s 

1.0 

51 

51 

52 

52 

53 

54 

56 

58 

62 

69 

83 

Table  6.11:  Number  of  iterations  for  the  two  region  scat¬ 
tering  ratio  sweep  using  SI,  DSA,  diamond 
difference,  S4  single-range  Gauss-Legendre 
angular  quadrature  in  XY  geometry. 


Left  Region 

Scattering  Ratio 

0.0 

0.1 

0.2 

:  0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

.2 

0.0 

2 

6 

8 

9 

9 

10 

10 

11 

11 

12 

13 

cd 

0.1 

6 

7 

8 

8 

9 

9 

10 

11 

12 

12 

13 

bO 

0.2 

7 

7 

7 

8 

9 

9 

10 

11 

12 

12 

13 

Th 

CD 

0.3 

8 

8 

8 

8 

9 

9 

10 

11 

12 

12 

13 

0.4 

8 

8 

8 

8 

9 

9 

10 

11 

12 

13 

13 

CD 

CO 

0.5 

9 

8 

8 

9 

9 

9 

10 

11 

12 

13 

14 

£ 

O 

0.6 

9 

9 

9 

9 

9 

9 

10 

11 

12 

13 

14 

‘5o 

CD 

0.7 

10 

10 

9 

9 

9 

9 

10 

11 

12 

13 

14 

0.8 

10 

10 

10 

10 

10 

10 

10 

11 

12 

13 

14 

2 

bp 

0.9 

11 

11 

11 

11 

11 

11 

10 

11 

12 

13 

14 

2 

1.0 

12 

12 

12 

12 

12 

12 

11 

11 

12 

13 

14 
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Table  6.12:  Number  of  iterations  required  for  conver¬ 
gence  in  the  time  independent,  XYZ  geom¬ 
etry  tests. 


Test  Problem 

DD 

DI 

Step 

Source  Iteration 

Unaccelerated  DSA 

Pure  Absorber 

2 

2 

2  2 

Pure  Scatterer 

37 

8 

410  11 

Infinite  Medium 

NC 

14 

421  388 

Constant  Source,  Single-Region 

11 

8 

17  8 

thus  the  total  number  of  iterations  in  the  “Three  Group  Absorber,  Upscatter”  test 
was  greater  than  what  PARTISN  required  (30  for  DI  and  18  for  unaccelerated  source 
iteration). 

The  multigroup  problems  were  tested  using  three  different  initialization  strate¬ 
gies  for  the  inward  angular  current  distribution  for  each  energy  group:  Isotropic  at 
every  outer  iteration  (Isotropic);  isotropic  for  the  first  outer  iteration  and  the  previous 
distribution  for  subsequent  outer  iterations  (Previous);  and  a  random  initialization 
for  every  outer  iteration  (Random).  As  shown  in  Table  6.14,  for  downscatter  prob¬ 
lems  there  is  no  advantage  between  Isotropic  and  the  Previous  initialization  methods. 
For  problems  that  have  upscatter,  the  Previous  method  performs  better  than  the 
Isotropic  method.  The  Random  initialization  method  demonstrates  that  convergence 
is  not  predicated  on  the  initial  distribution  when  solving  multigroup  problems. 


Table  6.13:  Total  number  of  inner  iterations  required  for 
convergence  in  the  time  independent,  multi¬ 
group,  slab  geometry  tests. 


Source  Iteration 

Test  Problem 

DI 

Unaccelerated 

DSA 

TSA 

Three  Group  Absorber,  Downscatter 

6 

6 

6 

6 

Three  Group  Alternating 

11 

135 

13 

34 

Three  Group  Absorber,  Upscatter 

21 

21 

9 

21 
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Table  6.14:  Total  number  of  inner  iterations  required  for 
convergence  in  the  time  independent,  multi¬ 
group,  slab  geometry  inward  angular  partial 
current  initialization  method  tests. 


Test  Problem 

Isotropic 

Previous 

Random 

Three  Group  Absorber,  Downscatter 

6 

6 

6 

Three  Group  Alternating 

11 

11 

13 

Three  Group  Absorber,  Upscatter 

24 

21 

24 

6.3  Time  Dependent  Verification  Tests 

The  agreement  between  PARTISN  and  DI  was  sufficient  to  demonstrate  that 
DI  can  be  used  as  the  transport  code  in  a  time  dependent  problem.  The  difference 
that  was  observed  in  the  two  source  ramp  tests  was  due  to  the  differences  in  the 
time  steps  and  where  the  time-varying  source  was  sampled.  The  number  of  iterations 
required  for  each  time  step  was  consistent  with  the  results  seen  in  the  time  independent 
testing,  e.g.  only  two  iterations  are  needed  when  the  synthetic  scattering  ratio  was 
small.  For  the  pure  scatter  case,  when  vAt  was  increased,  thereby  allowing  c  — »  c, 
the  number  of  iterations  increased.  Also,  when  the  synthetic  cross  section  was  large, 
the  number  of  iterations  was  generally  insensitive  to  changes  in  the  actual  total  cross 
section  because  the  synthetic  scattering  ratio  was  relatively  unaffected.  Both  DI  and 
PARTISN  exhibited  these  behaviors,  which  supports  the  conclusion  that  DI  performs 
in  an  equivalent  fashion  as  PARTISN. 
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VII.  Conclusion 


The  Distribution  Iteration  method  has  consistently  been  demonstrated  as  a 
reliable  replacement  for  source  iteration.  Prior  to  my  research,  there  had  been  no 
formal  verification  plan  for  distribution  iteration.  The  previous  work  by  Wager  and 
Prins  only  implemented  tests  designed  to  aid  in  the  development  of  the  algorithm. 
The  verification  plan  that  was  implemented  did  not  reveal  any  unexpected  difference 
between  the  results  produced  by  DI  and  accelerated  source  iteration. 

7. 1  Contributions 

In  my  research,  I  did  demonstrate  that  DI  can  be  used  to  solve  the  time  depen¬ 
dent  transport  equation  using  the  same  method  used  by  production  source  iteration 
codes.  As  demonstrated  in  this  dissertation,  there  is  no  optimal  combination  of  pa¬ 
rameters,  e.g.  time  step  and  cell  size,  for  producing  reliable  results.  The  choice  of  a 
time  step  becomes  more  complicated  with  multigroup  problems.  While  it  is  possible  to 
pick  a  time  step  and  still  obtain  a  reasonable  cell  optical  thickness  for  monoenergetic 
problems,  a  multigroup  problem  requires  a  transport  algorithm  that  tolerates  opti¬ 
cally  thick  cells  for  some  energy  groups.  My  research  demonstrated  that  DI  was  able 
to  produce  reasonable  results  for  optically  thick  time  dependent  problems  without  a 
fix- up  mechanism. 

One  of  the  design  goals  for  the  software  that  I  developed  as  part  of  my  re¬ 
search  was  to  be  able  to  support  new  research.  I  implemented  a  flexible  and  modular 
code,  which  also  incorporated  real  time  visualization.  Given  the  uniqueness  of  the  DI 
algorithm,  I  had  to  develop  the  visualization  methods. 

In  order  to  be  able  to  use  the  TIEL  benchmark  developed  by  Ganapol,  Mathews 
and  I  developed  the  matrix  albedo  for  time  independent  problems.  This  new  method 
provides  an  accurate  solution  for  the  effect  that  an  infinite  medium  will  have  on  the 
angular  flux  in  a  region  of  interest.  While  some  aspects  of  the  implementation  were 
straightforward,  I  implemented  all  boundary  conditions,  e.g.  specular  reflection,  in 
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matrix  form.  This  new  approach  provides  a  flexible  method  for  handling  difficult 
boundary  conditions,  such  as  asymmetric  angular  quadratures. 

Several  improvements  to  the  D1  algorithm  are  original  to  my  research.  The 
first  improvement  was  the  methodical  implementation  of  the  three  geometries.  This 
implementation  is  more  compact  and  maintainable  than  previous  versions  and  was 
instrumental  in  keeping  verification  plan  simple.  The  improvement  to  the  angular 
current  distribution  refinement  algorithm  converges  more  rapidly  than  the  previous 
algorithm.  The  new  algorithm  is  also  more  parallel  efficient  than  the  previous  version. 
Finally,  the  partial  current  solver  can  easily  be  adapted  to  different  linear  algebra 
solvers  and  has  improved  parallel  efficiency. 

7.2  Conclusions 

Wager  and  Prins  performed  initial  assessments  of  the  relative  performance  be¬ 
tween  DI  and  SI.  My  research  examined  a  broader  range  of  the  parameter  space  using 
both  DSA  and  TSA.  Based  on  my  research,  DI  consistently  performed  as  well  or 
better  than  either  DSA  or  TSA  in  slab  geometry.  In  XY  and  XYZ  geometries,  DI 
continues  to  exhibit  superior  performance  when  using  positive  spatial  quadratures. 
The  performance  of  DI  relative  to  DSA  and  TSA  supports  the  conclusion  that  DI  is 
a  powerful  replacement  for  source  iteration.  For  the  classes  of  problems  where  DSA 
loses  effectiveness  and  TSA  fails,  DI  exhibits  robust  performance  and  converges  in 
fewer  iterations. 

In  XYZ  geometry,  DI  does  demonstrate  a  pronounced  weakness  when  using  a 
spatial  quadrature  that  generates  negative  fluxes.  The  poor  performance  exhibited 
with  diamond  difference  is  strong  evidence  that  DI  should  be  used  with  a  positive  spa¬ 
tial  quadrature  such  as  SC  or  with  a  spatial  quadrature  that  is  negative  less  frequently 
than  DD,  e.g.  LC  or  LD.  I  did  implement  a  new  negative  angular  partial  current  dis¬ 
tribution  fix-up,  which  is  different  then  the  fix-up  that  PARTISN  performs,  however, 
it  was  ineffective  in  highly  scattering  problems  (c  >  0.7). 
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For  time  dependent  transport,  PARTISN  implements  an  adaptive  time  step 
algorithm  and  a  negative  flux  fixup.  In  addition,  PARTISN  enables  the  negative  flux 
fixup  algorithm  in  the  transport  equation.  Reverse  engineering  these  three  algorithms 
was  impractical  and  of  little  value.  The  comparisons  were  adequate  to  confirm  that 
my  time-dependent  DI  code  has  no  errors  worse  than  PARTISN’s. 

The  poor  performance  that  DI  exhibited  when  using  DD  is  not  a  failure  of  the 
DI  algorithm.  Mathews  [26]  demonstrated  the  unphysical  ray  effect  that  DD  exhibits, 
which  calls  into  question  the  reliability  of  the  results  generated  by  DD  when  negative 
fluxes  occur,  independent  of  the  algorithm.  A  key  advantage  that  DI  has  is  that 
new  spatial  quadratures  can  be  implemented  easily  while  implementing  a  new  spatial 
quadrature  in  DSA  is  not  trivial. 

Currently,  there  is  no  acceleration  implemented  in  the  DI  algorithm.  In  slab 
and  XY  geometry  there  was  little  motivation  to  add  any  acceleration  when  using  DI 
because  it  outperforms  DSA  and  TSA.  In  XYZ  geometry  my  research  shows  that 
acceleration  will  have  little  benefit  when  using  a  positive  spatial  quadrature.  There 
may  be  some  classes  of  problems  where  DI  may  benefit  from  a  convergence  accelerator. 
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Appendix  A.  Other  Transport  Codes 


A.l  TIMEX 

The  within-group  BTE  used  by  the  1972  TIMEX  code  [30]  is  originally  presented 
as  (N.B.  some  subscripts  have  been  eliminated  for  clarity) 


1  ( Tt1  -  Tl 

v  1  At 


Ai+ 1/2^+172  _  Ai-1/2^1+1/ 


WmVi 


+  =  SJmti  +  qm,i.  (A.l) 


TIMEX  1972  uses  a  cell-centered  indexing  scheme,  where  j  is  the  time  index,  i  is  the 
cell  index  and  m  is  the  ordinate  index.  Cell  faces,  therefore,  are  indicated  by  a  half 
step,  e.g.  i  +  1/2  is  the  cell  face  between  cells  i  and  i  +  1.  The  area  of  a  cell  face  is 
represented  by  A  and  the  cell  volume  by  V.  The  ordinates  and  weights  are  represented 
by  /i  and  w,  respectively.  The  variable  a  represents  the  curvature  coefficient.  The 
angular  flux  is  represented  by  ip,  S  is  the  scattering  and  fission  sources  from  the  j 
time  interval,  and  q  is  the  inhomogeneous  source.  Finally,  v  is  the  group  velocity  and 
At  is  the  time  step. 

The  combination  of  A,  V  and  a  allows  the  TIMEX  1972  code  to  solve  three 
different  one-dimensional  geometries-slab,  cylindrical,  and  spherical.  The  curvature 
coefficient  is  defined  recursively  as 


^m+ 1/2  1/2  pr 


Ai~  1/2) , 


(A.2) 


using  ay/2  =  a^+ 1/2  =  0,  where  N  is  the  number  of  ordinates,  as  starting  conditions 
(N.B.  the  ordinate  index,  m  starts  at  1).  The  definitions  for  A  and  V  are  shown  in 
table  A.l.  Note  that  a  recursive  definition  is  used  for  the  area  elements  in  spherical 
geometry.  This  was  done  to  improve  accuracy  for  cells  near  the  center  of  the  sphere. 
Also,  in  cylindrical  geometry  the  ordinates  represented  by  pm  and  wm  is  actually  a 
two-dimensional  angular  quadrature  of  /j  and  £. 
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Table  A.l:  Definition  of  geometric  functions  used  by 
TIMEX  1972. 


Geometry 

Variable 

Cell  Face  Area  (A) 

Cell  Volume  (V) 

Slab 

Cylindrical 

Spherical 

25+1/2 

r t+l/2 

r  *+1/2 

1 

27rr  j+i/2 

2D  /l  . 

25+1/2  —  21j_  1/2 
^  (c+1/2  -  ri~l/2) 

47T  /  3  ,r3  \ 

3  y  *+1/2  '  i— 1/2  y 

^i+1/2  ri  —  l/2  1  1 /2 

TIMEX  (1972)  uses  a  cell  average  relationship  between  the  cell-centered  flux 
and  the  edges  in  both  space  and  angle 


P+1  - 

Ah + Cb 

(A.3) 

i,m 

2 

P+1  _ 

A+V + cc 

(A.4) 

i,m 

2 

These  relationships  are  used  eliminate  the  i  +  1/2  and  m  +  1/2  for  fim  >  0  directions 
and  i  —  1/2  and  m  —  1/2  for  /im  <  0  directions  and  obtain 

i^vAt  +  ^m^i+1/2  - w  ~ ^  a*^^  ^i+1  =  (Ai+ 1/2  +  Ai- 1/2)  V£u/2 

1/2 )  ^+-i/2 + (^)  ^ + + (A-5) 

The  recursion  relationship  for  the  curvature  coefficient  is  used  to  obtain 


Vi  +  (Aw/2  +  +_1/2)  +  a’"+1/2.+  Q""1/2  +  *V5 )  A' 


uAf 


wr, 


—  Hm  {Ai+ij 2  +  A-!/2)  1/2  + 


<*m+l/2  +  am-l/2  \  ,/j+l 

Vm- 1/2 


Wr 


uAf 


$  +  (A. 6) 
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Because  my  research  is  focused  on  Cartesian  geometry,  the  BTE  simplifies  to  (substi¬ 
tuting  A  Xi  =  Xi+ 1/2  —  Xi- 1/2  and  dividing  by  A  ay) 


+  2 


vAt  Ax . 


+  * )  d+1  = 


A+1  i  1 

"1//2  uAf 


+  Sln,i  +  % 


(A.7) 


Substituting  the  relationship  for  the  midpoint  spatial  average  and  algebraically  rear¬ 
ranging  the  terms  yields 


vAt 


+  <?i  )  ^  +  d 


wd+i 

V/+1/2 


-C1 


1/2 


Aay 


uAt 


+  S'm.i  +  ft 


m,2  • 


(A.8) 


The  use  of  the  diamond- difference  spatial  discretization  can  result  in  a  negative 
outgoing  angular  flux.  Using  the  midpoint  spatial  average  to  express  + 1  in  terms 
°f  ^i+1/2  and  ^i—i/2  yiclds 


f  1  ,  \  ^i+l/2  +  ^S/2  ,  ^i+1/2  -  C/ 2 

+  aV  - 2 - +  ^ - Aay - 


vAt 


H  +  ■Sm,  +  ft 


m,i • 


(A.9) 


Solving  for  i/fff  2  in  terms  of  V;ff/2  yields  the  equation 


J+1 


w^'+i  — 

Vi+l/2  — 


i  Air 

1  2/x 


A  +  di) 


^i-l/2  + 


1  +  If  (<u  +  ds) 


1/2 


1  +  if  (a*  +  si*) 


(A.10) 


The  /2  coefficient  will  be  negative  when 


(JiAx  + 


Ax 

vAt 


>2|/i|. 


(A.ll) 


When  <t i  +  -A-  >>  1,  the  source  term  will  be  negligible  and  either  a.iAx  or  ^  (or 
both)  will  be  large.  Consider  the  case  where  At  — >  0, 


3+ 1 
i+1/2 


/  i  /  7+1 

aW  ~  A 


1/2- 


(A.12) 


106 


Thus,  ^\X\/2  has  the  potential  for  being  negative  if  the  cell  centered  flux  is  small, 
particularly  when  /i  is  small. 

Removing  the  discretization  of  the  spatial  derivative  from  equation  A. 8  yields 


1  d  \ 

^At+a‘  +  ^mdx) 


dd  =  +  sLi  +  (i 


m,f 


(A. 13) 


Similarly,  the  discretization  of  the  time  derivative  can  be  removed,  which  yields 


fid  d\ 

{vdt  +  a‘  +  ,lmai) 


t 


_  gj 


+  Qn 


(A. 14) 


In  operator  notation,  this  is  equivalent  to 


+  l¥+1  =  S^>j  +  q, 

ot 


(A. 15) 


where 


d 


(A. 16) 


A. 2  PARTISN 

PARTISN  is  based  on  the  implicit  in  time  form  of  the  BTE.  Starting  from  the 
equation  (Alcouffe  [4]) 

(l§t  +  il'v  +  a)i’  =  s  (A-17) 

and  using  a  central  difference  for  the  time  derivative  yields 

—  (Vd+1/2  -  ^”1/2)  +  Cl  ■  +  a^j  =  Sj.  (A. 18) 

vAt  v 


The  above  equation,  which  is  equivalent  to  the  explicit  midpoint  method,  has  second 
order  local  truncation  error  in  time,  i.e.  0((Af)2).  The  equation  can  be  returned  to 
an  implicit  form  by  eliminating  ?/d +1 /2.  Consider  the  midpoint  time  average 


i\)i 


^■+1/2  +  ^1-1/2 
2 


(A. 19) 
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The  midpoint  time  average  can  be  rewritten  to  solve  ?/d+1/2  as 


ipj+1/2  =  2^  -  ip j~lf 2.  (A. 20) 

We  can  eliminate  ipi +1 /2  by  using  the  midpoint  time  average  obtaining, 

-4-  (^'  -  ^'-1/2)  +  fi  •  v^'  +  =  Sj.  (A. 21) 

vAt 

The  above  equation  is  now  a  backwards  difference  and  has  first-order  local  truncation 
error.  Regrouping  the  terms  in  a  similar  fashion  to  the  TIMEX  1976  approach  yields 

O  •  Wipj  +  (a  +  ^  (A.22) 

The  flux  at  the  j  +  1/2  time  step  is  recovered  by  using  equation  (A. 20)  combined  with 
the  result  for  ?/d  obtained  from  equation  (A. 22). 

The  time  stepping  algorithm  implemented  by  PARTISN  is  shown  in  figure  A.l. 
Because  PARTISN  was  adapted  from  the  DANTSYS  transport  code,  using  half  time 
steps  for  the  extrapolated  flux  was  notationally  convenient.  The  half  time  steps  are 
computed  by  the  midpoint  extrapolation  and  the  integer  intervals  are  computed  by 
the  transport  algorithm.  The  notation  is,  however,  awkward  and  can  be  expressed 
instead  as 

12  •  W  +  (v  +  V  =  Sj  +  (A. 23) 

where  At'  =  At/2.  The  extrapolation  step  then  becomes 

ipj+1  =  2 ipj  -  (A. 24) 

It  is  important  to  note  that  the  truncation  error  associated  with  the  i/d  angular  flux  is 
first  order  because  it  is  backward  Euler-there  is  no  apparent  benefit  in  the  truncation 
error  to  using  this  approach.  The  improvement  in  the  truncation  error  is  in  the 
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extrapolation  step.  Doing  a  Taylor  series  expansion  on  ?/d"fl  and  rjp  1  yields 


^j+1  =  +  At' ^ 

at 

vj  1  =  i)j  -  At'0'' 
at 


(At')2  d2ip 
+  2  dt2 

(At')2  d2^ 
+  2  <9f2 


+  0((Af')3) 
-  0((Af')3). 


(A. 25) 
(A. 26) 


Substituting  in  the  Taylor  series  expansion  for  ?/T  1  into  the  extrapolation  step  yields 


?/d+1  = 


^  +  At 


dt  j 


(At')2  d2i> 
~^~~dt2 


+  0((At')3). 


(A.27) 


Note  that  the  extrapolation  step  results  in  a  ?/d+1  that  is  accurate  to  (At')2  term. 
Thus,  this  method  will  yield  second-order  local  truncation  error  only  with  the  extrap¬ 
olated  fluxes  and  not  with  the  midpoint  fluxes. 

PARTISN,  however,  will  only  output  fluxes  (scalar  or  angular)  from  the  mid¬ 
point  time  and  not  the  extrapolated  end  of  time  step  flux.  Thus,  the  results  from 
PARTISN  are  only  accurate  to  first  order  even  though  the  algorithm  is  accurate  to 
second  order  internally.  PARTISN  solves  (A. 22),  which  is  effectively  Backwards  Euler, 
in  the  transport  solver  and  then,  in  a  separate  state,  applies  the  extrapolation  step. 
One  advantage  to  this  approach  is  that  the  global  error  will  scale  by  (At)2  rather 
than  by  At  if  backward  Euler  was  used. 


Initialize  "01/2  to  the  starting  angular  flux 
for  j  =  1  to  Number  of  Time  Steps  do 
Compute  from  the  transport  equation 
Use  the  midpoint  average  to  extrapolate  to  ?/d+1/2 
Increment  the  time  index 

The  extrapolated  flux  becomes  the  new  '0J~1/2  angular  flux 

end  for 

Algorithm  A.l:  Semi-Implicit  Time  Stepping  Algo¬ 
rithm 

In  order  to  conduct  the  testing  to  support  my  research,  two  changes  to  the 
PARTISN  code  were  made.  The  Erst  change  was  to  have  PARTISN  preserve  the 
scalar  flux  (RTFLUX)  for  each  time  step.  By  default,  PARTISN  discards  the  scalar 
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flux  from  each  time  step  and  there  is  no  parameter  hie  option  to  change  that  behavior. 
The  second  change  was  to  correct  an  erroneous  time  stamp  in  the  time  dependent 
scalar  flux.  PARTISN  was  writing  the  end  of  time  step  rather  than  the  midpoint 
time,  which  was  the  scalar  flux  that  was  being  written  to  the  hie. 


110 


Appendix  B.  Transport  Coefficients 


Wager  [32]  and  Prins  [28]  derived  the  transport  coefficients  in  slab  and  XY 
geometries  for  various  spatial  quadratures.  In  the  following  sections,  I  will  present 
the  derivation  of  the  Step  method  spatial  quadrature  in  XY  and  XYZ  geometries  and 
Diamond-Difference  in  XYZ  geometry. 

The  derivations  are  shown  using  the  time  independent  form  of  the  BTE.  If 
time  dependent  coefficients  are  needed  for  the  synthetic  time  form  of  the  BTE,  both 
the  cross  section  term,  a  and  the  average  source  term,  Sa,  are  replaced  with  the 
corresponding  synthetic  form. 


B.l  Step  Method 

The  Step  method  is  defined  by  the  assumption  that  the  angular  flux  is  constant 
within  a  cell  (supressing  the  cell  index  number  for  clarity),  specifically 


and 


H  >  0 

(B.l) 

n  <  0 

(B.2) 

IpT  =  -0A  ?7  >  0 

(B.3) 

Tb  =  V  <  0, 

(B.4) 

where  /i  and  rj  are  the  direction  cosines  from  the  angular  quadrature.  This  spatial 
quadrature  is  very  inaccurate,  but  it  does  have  the  benefit  of  being  positive  and  useful 
when  debugging.  The  derivation  begins  with  the  zeroth-moment  balance  equation 


fpR  -  fpL 


+  V 


Ay 


+  0-fA 


Sa- 


(B.5) 
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Kf 

Kf 

K^s 


a 

€y  +  O'  +  1 

1 

6y  +  Oi  +  1 

Ay 

- - - -  and 

£y  +  o;  +  1 
Ay 
V 

£y  +  +  1 


(B.12) 

(B-13) 

(B-14) 

(B-15) 
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The  coefficients  for  the  contributions  to  the  X  faces  (left  and  right  faces)  are 


k2  = 
K  = 

k°e- 

k?e 


a 


£v  +  OL  +  1 


£y  +  OL  +  1 
Ay 

and 


£y  +  OL  +  1 
Ay 
V 

£v  +  OL  +  1 


(B.16) 

(B-17) 

(B-18) 

(B.19) 


The  coefficients  for  the  contributions  to  the  Y  faces  (bottom  and  top  faces)  are 


KS  = 


Kv  = 


a 


= 


£v  +  OL  +  1 


+  a  +  1 

Ay 

and 


K°E  = 


ev  +  o;  +  1 

Ay 
>1 

£v  +  o;  +  1 


(B.20) 

(B.21) 

(B.22) 

(B.23) 


The  derivation  for  the  X-Y-Z  geometry  is  similar,  with  the  solution  for  the  cell 
average  angular  flux  as 


SA^f  +  Oii>L  +  'lpB  +  Pilp 
£y  +  OL  +  (3  +  1 


(B.24) 
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where  (3  =  gv/gz  and  ez  =  crAz /£.  The  transport  coefficients  for  the  contributions  to 
the  cell  average  angular  flux  are 


Kf  = 

a 

(B.25) 

Gy  +  OL  +  f3  +  1 

Kf  = 

1 

(B.26) 

Gy  +  OL  +  f3  +  1 

Kf  = 

p 

(B.27) 

Gy  +  OL  +  f3  +  1 

Ay 

K^s  = 

Tnrl 

(B.28) 

c  1 1 1  v  1 

€y  +  OL  +  (3  +  1 

Ay 

K^e  = 

V 

(B.29) 

Gy  +  01  +  (3  +  1 

The  coefficients  for  the  contributions  to  the  X  faces  (left  and  right  faces)  are 


k2  = 

a 

(B.30) 

Gy  +  OL  +  f3  +  1 

kS  = 

1 

(B.31) 

Gy  +  OL  +  f3  +  1 

k2  = 

p 

(B.32) 

Gy  +  Of  +  (3  +  1 

Ay 

Kf  = 

>1  nnrl 

(B.33) 

dllU. 

+  cr  +  (3  +  1 

Ay 

Kf  = 

V 

(B.34) 

Gy  +  OL  +  (3  +  1 
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The  coefficients  for  the  contributions  to  the  Y  faces  (bottom  and  top  faces)  are 


= 

a 

(B.35) 

Gy  +  OL  +  f3  +  1 

Kv  = 

1 

(B.36) 

Gy  +  OL  +  f3  +  1 

p 

(B.37) 

Gy  +  OL  +  (3  +  1 

Ay 

K™  = 

v  nnrl 

(B.38) 

dllU. 

Gy  +  CY  +  (j  +  1 

Ay 

K°E  = 

V 

(B.39) 

Gy  +  Ol  +  (3  +  1 

The  coefficients  for  the  contributions  to  the  Z  faces  (back  and  front  faces)  are 


K®1 = 

a 

(B.40) 

Gy  +  OL  +  f3  +  1 

KS  = 

1 

(B.41) 

Gy  +  OL  +  f3  +  1 

k2  = 

p 

(B.42) 

Gy  +  OL  +  (3  +  1 

Ay 

v  nnrl 

(B.43) 

el!  1  v  1 

Gy  +  a  +  p  +  l 

Ay 

k°e  = 

V 

(B.44) 

Gy  +  Ot  +  (3  +  1 

B.  2  Diamond  Difference 

The  Diamond  Difference  spatial  quadrature  applies  the  assumption  that  the 
angular  fluxes  are  related  to  the  cell  average  angular  flux  by 

tpL  +  tpR  ~  IpB  +  lpT  —  fpF  +  tpP  —  (B.45) 


as  an  auxiliary  equation  to  the  cell  balance  equation 


tpR  ~ 


tpT-tpB  ,  F  ~  fpP 
+  V — y—  +  f - 


Ay 


Az 


+  (Tip  a 


SA. 


(B.46) 
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Using  Equation  B.45  to  solve  for  ipR,  ipp,  and  f  in  terms  of  ipL,  ipB,  ipP  and  ip  a 
yields 


ipR  =  2  ip A  -  ipL , 
ipT  =  2  ip  a  —  ipB  and 
ipF  =  2 ip A  -  ipp. 


(B.47) 

(B.48) 

(B.49) 


Using  the  above  equations  to  eliminate  ipB,  ipp  and  ipP  and  grouping  like  terms  yields 


a + +  + 2~k)  *A  ~  Sa + 2£v’1 + 2£/b + 2h'pp'  (B-50) 


Let 


ey 


a  Ax 

1 

U 

aAy 

V 

crAz 


ez  = 


£ 


a  =  —  and 


y 

ez 


(B.51) 

(B.52) 

(B.53) 

(B.54) 

(B.55) 


The  ratios  a  and  /3  was  chosen  to  be  consistent  with  the  X-Y  diamond  difference 
equation  derived  by  Mathews  [26].  Multipling  Equation  B.50  by  A y/rj  yields 


Av 

(cy  +  2<y,  +  2  +  2/3)  ip  a  —  Sa - h  2  oiipp  +  2  ipB  +  2  phpp. 


Solving  for  the  cell  average  scalar  flux  provides 


(B.56) 


SA^  +  2  aiPL  +  2ipB  +  2phpP 
r  ey  +  2a  +  2  +  2/3 


(B.57) 
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Solving  for  ipR,  ipT  and  tfjP  yields 


'tpR 


^ T 


'tpF 


2  SA^  +  (2a  -ey-2-  2  /3)^l  +  +  4/3^ 

6y  +  2a  +  2  +  2/3 

2 SA^  +  4a^L  +  (2  -  ey  -  2a  -  2/3)^  +  4^P 
ey  +  2a  +  2  +  2/3 

2^^  +  4a^L  +  4 tl>B  +  (2/3  -ey-2a-  2)^P 
ey  +  2a  +  2  +  2/3 


(B.58) 

(B.59) 

(B.60) 


Using  the  above  solutions,  the  transport  coefficients  for  the  contributions  to  the  cell 
average  angular  flux  are 


Kf  = 

2a 

(B.61) 

Gy  +  2a  +  2/3  +  2 

Kf  = 

2 

(B.62) 

ey  +  2a  +  2/3  +  2 

y 

Kf  = 

2/3 

(B.63) 

ey  +  2a  +  2/3  +  2 

Ay 

K*s  = 

v  nnd 

(B.64) 

Cllivl 

+  2a  +  2/3  +  2 

Ay 

KyE  = 

V 

(B.65) 

(B.66) 

Gy  +  2a  +  2/3  +  2 
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The  coefficients  for  the  contributions  to  the  X  faces  (left  and  right  faces)  are 


= 

2a- gv-2-2  (3 

(B.67) 

XX 

Cy  +  2  a  +  2/3  +  2 

= 

4 

(B.68) 

xy 

Gy  +  2a  +  2(3  +  2 

k£  = 

4  (3 

Gy  +  2a  +  2(3  +  2 

2^i 

(B.69) 

K°s  = 

- ^ ^  and 

(B.70) 

X 

Gy  +  2a +  2(3 +  2 

2Ay 

k°e  = 

V 

Gy  +  2a  +  2(3  +  2 

(B.71) 

(B.72) 

The  coefficients  for  the  contributions  to  the  Y  faces  (bottom  and  top  faces)  are 


= 

4  a 

(B.73) 

yx 

Gy  +  2a  +  2(3  +  2 

= 

2  -  Gy  -  2a  -  2(3 

(B.74) 

yy 

Gy  +  2a  +  2(3  +  2 

= 

A(3 

(B.75) 

yz 

Gy  +  2a  +  2(3  +  2 

2+l 

K°s  = 

- o — — o  and 

(B.76) 

y 

Gy  +  2a  +  2(3  +  2 

2+y 

K°E  = 

y 

(B.77) 

y 

Gy  +  2a +  2(3 +  2 

(B.78) 
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The  coefficients  for  the  contributions  to  the  Z  faces  (back  and  front  faces)  are 


= 

4  a 

(B.79) 

6y  +  2a  +  2/3  +  2 

= 

4 

(B.80) 

zy 

6y  +  2  a  +  2/3  +  2 

= 

2/3  —  ey  —  2a  —  2 

(B.81) 

zz 

6y  +  2a  +  2/3  +  2 

2As 

Kos  = 

= - - — — -  and 

(B.82) 

z 

ey  +  2a  +  2/3  +  2 

2Ay 

koe  = 

y 

(B.83) 

z 

£y  +  2a  +  2/3  +  2 

(B.84) 
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Appendix  C.  Boundary  Currents 


Currents  that  are  incident  on  the  boundary  of  the  region  of  interest  can  be 
classified  by  the  treatment  of  the  source  term.  Consider  the  following  simplified, 
linearized  BTE  for  a  pure  absorber  (c  =  0)  medium, 

dib 

+  cnl)  =  Q(x,y).  (C.l) 


For  illustrative  purposes,  we  can  further  simplify  this  example  by  assuming  that  the 
incident  current  is  from  a  source  on  the  left  boundary  and  there  is  no  source  of 
neutrons  internal  to  the  problem,  thus  all  the  neutrons  are  moving  in  the  y  >  0 
direction.  A  current  source  on  the  left  can  be  expressed  as 


Q(x,y) 


{Q(y)  when  x  =  0 
0  otherwise. 


(C.2) 


We  can  group  the  current  sources  into  two  categories  based  on  whether  or  not  there 
is  any  y  dependence  in  the  source. 

A  Lambertian  or  diffuse  source  [6]  is  one  that  has  the  same  angular  flux  when 
viewed  from  any  angle,  e.g.  the  detectors  at  A  and  B  as  shown  in  Figure  C.l  both 
observe  the  same  angular  flux  from  the  differential  surface  element  dA  Thus,  a  Lam¬ 
bertian  source  has  a  uniform  flux  per  steradian-Q(/i)  =  constant.  Solving  equation 
C.l)  for  positive  y  yields 

y)  =  Qe~xa/,i.  (C.3) 

Integrating  over  all  y  yields  the  scalar  flux,  0, 

=  j  (e~xa  -  xaT(0,  xcr))  ,  (C.4) 


where  T (a,  z)  is  the  upper  incomplete  Gamma  function  [1]  defined  as 


r(a,  z) 


(C.5) 
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A 


As  the  left  boundary  is  approached,  <f)(x)  — >  Q/2  .  In  terms  of  current,  j,  along 
an  ordinate  a  Lambertian  source  will  have  the  familiar  cosine  dependence:  j  =  pQ. 
Integrating  over  p  >  0,  the  rightward  (positive)  partial  current  on  the  left  face  is  Q/4 
(using  the  normalization  scheme  introduced  by  Lewis  and  Miller  [22]). 

There  is  one  specific  type  of  source  of  the  second  type  that  is  useful  in  bench¬ 
marking  transport  codes-those  that  produce  a  uniform  current  source,  viz.  j  = 
constant.  This  source  is  characterized  by  the  uniform  current  flowing  from  the  sur¬ 
face,  hence  the  appellation  isotropic  surface  source  (ISS).  The  angular  flux  for  this 
type  of  source  is 

(C.6) 

p 

with  the  scalar  flux  (for  positive  x)  being 

<f>(x)  =  ®r(0,xa).  (C.7) 

As  the  left  boundary  is  approached,  <f>{x)  — >  oo.  The  rightward  partial  current  on  the 
left  face  is  Q/2. 
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The  key  differences  to  note  between  the  Lambertian  and  the  isotropic  surface 
sources  are  the  partial  currents  (Q/ 4  versus  Qj 2,  respectively)  and  the  finite  versus 
infinite  scalar  flux  at  the  current  surface,  as  shown  in  Figure  C.2.  Based  on  these  two 
distinctions,  one  can  quickly  verify  how  a  particular  method  has  treated  an  incident 
source.  For  example,  PARTISN  utilizes  a  Lambertian  current  source. 


Figure  C.2:  Lambertian  Illumination  and  Isotropic  Sur¬ 
face  Source  Scalar  Fluxes 


In  the  case  of  Ganapol’s  TIEL  benchmark,  the  angular  flux  in  the  pure  absorber 


case  is 


t^H(q 


2/i 


(C.8) 


where  H(z )  is  the  Heaviside  step  function.  Integration  over  /i  yields  the  scalar  flux 


4>{x)  =  -T(0,xa)H(xa). 


(C.9) 


We  can  conclude  that  the  TIEL  benchmark  will  behave  identically  to  the  isotropic 
surface  source  near  the  origin.  Figure  C.2  also  shows  the  plot  of  the  scalar  fluxes  from 
the  TIEL  benchmark. 
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Glossary 


error  A  defect  due  to  an  incomplete  implementation  of 

the  model  and  is  not  due  to  lack  of  knowledge.,  61 

explicit  method  A  numerical  method  for  solving  differential  equa¬ 

tions  where  only  current  time  step  values  are  used 
to  compute  extrapolated  values,  14 

The  degree  to  which  a  model  is  an  accurate  repre¬ 
sentation  of  physical  process.,  60 
A  numerical  method  for  solving  partial  differential 
equations  where  all  derivatives  are  represented  with 
finite  differences,  14 

implicit  method  A  numerical  method  for  solving  differential  equa¬ 

tions  where  values  from  multiple  time  steps,  in¬ 
cluding  extrapolated  values,  are  used  to  compute 
extrapolated  values,  15 

isotropic  surface  source  A  current  source  where  the  current  along  each  or¬ 
dinate  is  constant.  See  Lambertian  source.,  120 

Lambertian  source  A  current  source  which  has  uniform  radiance  per 

steradian,  119 

positive  method  A  spatial  quadrature  that  will  unconditionally  gen¬ 

erate  positive  fluxes,  79 


fidelity 


fully  discretized 
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semi-discretization 


semi-implicit 


stiff  equation 


uncertainty 


validation 


verification 


A  numerical  method  for  solving  partial  differential 
equations  where  only  the  spatial  derivatives  are  rep¬ 
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